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Vector Operations 








This program provides solutions to the most common vector operations, 
such as addition, subtraction, dot and cross products, included angle, 
multiplication of a vector by a scalar, finding the length of a 
vector, and determining unit vectors. The program will allow vectors 
to be entered in either rectangular (x, y, z) or cylindrical (r, 
theta, z) coordinates, and will also display the result in either 
form. Conversion between these two formats is also an option. A useful 
feature of this program is the ability to chain vector operations by 
using the result from one operation as an operand for the next one. 


Let Vl = alX + b1Y + clZ and V2= a2X + b2Y + c2Z 


vl + V2 


(alta2)X + (b1l+b2)Y + (cl+c2)Z 
Vl - V2 = (al-a2)X + (bl-b2)Y + (cl-c2)2Z 
Vl X V2 = (bl*c2-b2*c1)X + (cl*a2-cl*c2)Y + (al*b2-a2*bl)Z 
Vl. V2 = (al*a2)X + (b1l*b2)Y + (cl*c2)2 
V1l*(s) = (al*s)X + (b1l*s)Y + (cl*s)Z 
norm(Vl) = (al*al + bl*bl + cl*cl)*.5 
unit vector(Vl) = (al/norm(Vl))X + (bl/norm(v1))Y¥ + (cl/norm(v1l))Z 


included angle 
between Vl & V2 


(norm(V1)*2 +norm(V2) *2-norm(V3) ~2) 


" 
» 
Q 
° 
n 
———1 


2*norm (V1) *norm(V2) 


V3 = vector between points defining Vl and v2 


page l 


User Instructions 


Comments Input Display 
[-----------------~------------- [----~--------- [--------------------- I 


1) Run the program. INPUT R/C ? 


2) The user may enter vectors 
in either rectangular or 
cylindrical form. Press [R] 
to enter x, y, z components 
of the vector or [C] to enter 
xr, theta and z components of 
the vector. If [C] is chosen, 
user flag 1 will be set as 
a visual reminder of the 
input format required. {R] or [C] DISPLAY R/C ? 


3) If the result of a vector 
operation is a vector, the 
program will display the 
result in either rectangular 
or cylindrical form. Press 
[R] to see the x, y, and z 
components of the result, or 
press [C] to see the result 
expressed in terms of r, 
theta and z. If [C] is 
chosen, user flag 2 will be 
set as a visual reminder of 
the display format. [R] or [C] A,S,X,1,D,N,M,U,C,F,Q? 


4) The program is prompting for 
the vector operation to be 
performed. The options are 
described below and grouped 
according to their input 
requirements. Vl represents 
the first vector entered, 
and V2 the second (if 
necessary). 


RRKKKKKKKKKRKKKKKKKKKRKKKRKKKRK KKK 


* SINGLE VECTOR OPERATIONS = 
KKK KKK RERKKEKEEEKRK KE KEKE 


N: calculate the norm or 
magnitude of a vector and 
return the scalar result. [N] 


M: multiply a vector by a 
scalar and return the 
vector result in the 
proper format. (M] 
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User Instructions 


Comments Input Display 


Us calculate the unit vector 
with the same direction 
as the input vector and 
with a norm of 1. Return 
the vector result in the 
proper format. [U] 


Cz: convert the input vector 
from rectangular to cyl- 
indrical or vice versa. 
This option requires that 
the input and display 
formats chosen by the user 
correspond to the desired 
conversion. If they do 
not, the message 'I/0 
INCORRECT’ will appear and 
the program will again ask 
for input and display 
formats. [C] 


For these single vector 
Operations, one of the 
following two sets of input 
Prompts will appear, depend- 
ing on the input mode 
selected. 


* rectangular * <x> [ENDLINE] y= 
* coordinates <y> [ENDLINE] z= 
<z> [ENDLINE] 


* 


* cylindrical * <r> [ENDLINE] theta= 
* coordinates * <theta>[ENDLN] z= 
<z> [ENDLINE] 


REKKKKKKKKKKKKKEKKKKKKKKRK KKK KKK 


* TWO VECTOR OPERATIONS = 
RHEKKKEKKKKKEKRKEKRKKKEKRKKKRKKKKERK RK RK 


A: add Vl and v2 and return 
the result in the proper 
format. [A] 
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User Instructions 
Comments Input Display 
J ------- ~~ - nr eee J -------------- J ------ - - ee ene I 
S: subtract V2 from Vl and 
return the result in the 
proper format. {s] 
X: calculate the cross 
product Vl x V2 and 
return the result in 
the proper format. {[X] 
I: calculate the included 
angle between V1 and V2 
and return the result 
in degrees. [I] 
D: calculate the dot or 
scalar product V1.V2 and 
return the scalar result. [D] 
For these two-vector opera- 
tions, one of the following 
two sets of input prompts 
will appear, depending on 
the input mode selected. 

VECTOR 1: x 
<xl>[ENDLINE] VECTOR 1: y= 
<yl>[ENDLINE] VECTOR l: z= 

* rectangular * <zl>[ENDLINE] VECTOR 2: x= 

* coordinates * <x2>[ENDLINE] VECTOR 2: y= 
<y2>[ENDLINE] VECTOR 2: z= 
<z2> [ENDLINE] 

VECTOR 1: r= 
<r1> [ENDLINE] VECTOR 1: theta= 
<thetal>[ENDLN] VECTOR 1: z= 

* cylindrical * <zl>[{ENDLINE] VECTOR 2: r= 
* coordinates * <r2>[ENDLINE] VECTOR 2: theta= 


KRkKkKRKKRKKKKKKRKKKKRKRKEKKRKRKKKRKKRKKEKEK 


* 


OTHER OPTIONS 


* 


kKkKkKRKKKKKKKRKKKKKKKKKKRKKKKKKKKKKEE 


F: 


allows the user to change 


input and/or display 
formats from cylindrical 
to rectangular or vice 
versa. 


<theta2> [ENDLN] VECTOR 


<z2>[ENDLINE] 


[F] 
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<go to step 2> 


User Instructions 


Comments Input Display 
I[--------------~--------~-~------ T---------+----- T-----------~--------- 1 


Q: exit program. [Q] > 


5) After an operation is chosen, 
the program will return the 
result in one of the formats 
listed below, depending on 
whether it is a vector or 
scalar. Pressing any key 
(except [ON]) will display 
the next component or move 
on to step 6. ( 


KKKKKKKKREREKKEKREREKKKKKE KKK KKK KKK 
mn VECTOR RESULTS 7” 


RkeKKKKKKKKKRKKKKKKRKKKKKKRKRKKKKKKKE 


[A],{S],(X),{M),{U),{[Cc) : 
These options will return a 
vector result in either 
rectangular or cylindrical 
coordinates, depending on 
the display option chosen 
in step 3. 


x= 
* rectangular * <any key> y= 
* coordinates * <any key> Z= 
<any key> <go to step 46> 
r= 
* cylindrical * <any key> theta= 
* coordinates * <any key> z= 
<any key> <go to step A> 
KRRKRKKKKKKKKKKKRKKKRKKKKKKR KKK KKK 
* SCALAR RESULTS * 
RARKKKKKKKKEKRKEKKKKKKKKRKR RK aK KKK 
{I] INCLUDED ANGLE = 
<any key> <go to step > 
[D] DOT PRODUCT = 


<go to step A> 


(N] NORM = 
<go to step 6> 
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User Instructions 


Comments Input Display 
[--------~----------------------- I-------------- I-------~-------------- I 


6) At this point, the user is 
given the option of exiting 
the program or running it 


again. RUN AGAIN (Y/N) ? 
To exit program, press 'N'. [N] > 
To run again, press 'Y'. [Y] <go to step 7> 


7) If the result of the last 
operation was a scalar, 
the program will return to 
step 4. 


If it was a vector, the 

program provides the option 

of using it in subsequent 

calculations. USE RESULT (Y/N)? 


To discard the result, press 
'N'. To use the result, 
press 'y'. (N] or [Y] <go to step 4> 


If 'Y' is chosen and the 
next operation requires 

one vector, the program will 
skip the prompt for the 
components and display the 
result. If two vectors are 
required, the program will 
assume the result is the 
first vector and will ask 
only for the second vector. 


\ 


REFERENCES: 
Salas, S. and Hille, E.; "Calculus", Xerox. 


Hudson, R.; “The Engineer's Manual", Wiley and Sons, Inc. 
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EXAMPLES 


A) Find the cross product of (3,5,8) and (4,8,1), and then calculate 
the angle the result vector makes with the x-axis (1,0,8). 


B) Find the length (norm) of the vector (12.66,-4.5,-7). 


C) Convert (-4,7,%) to cylindrical coordinates. 


Comments Input Display 
[------------------------------- [---~----------- [--------------------- I 
1) Run the program. [RUN] INPUT R/C ? 

2) Specify input in rectang- 
ular coordinates. {R] DISPLAY R/C ? 
3) Specify output in rectang- 
ular coordinates. {R] A,S,X,1,D,N,M,U,C,F,Q? 
4) Select cross product option. [X] VECTOR 1: x= 
5) Enter x, y, z for vectors. 3 [ENDLINE] VECTOR 1: y= 
5 [ENDLINE] VECTOR 1: z= 
8 [ENDLINE} VECTOR 2: x= 
4 [ENDLINE] VECTOR 2: y= 
8 [ENDLINE] VECTOR 2: z= 
L{ENDLINE] 
6) Display result. x=5 
<any key> =-3 
<any key> z=-26 
7) Run program again. <any key> RUN AGAIN (Y/N)? 
8) Use result. {Y] USE RESULT (Y/N)? 
9) Calculate included angle. {Y] A,S,X,1,D,N,M,U,C,F,Q? 
10) Enter second vector. [I] VECTOR 2: x= 
1 [ENDLINE] VECTOR 2: y= 
@ [ENDLINE] VECTOR 2: z= 
@ [ENDLINE] 
11) Display result. ANGLE = 76.1130063355 
12) Run program again. <any key> RUN AGAIN? (Y/N) 
13) Calculate norm. {YJ A,S,X,1,D,N,M,U,C,F,Q? 
14) Input vector. [N] x= 


12.66{ENDLINE] y= 
-4.5 [ENDLINE] z= 
-7 {ENDLINE] 

15) Display result. NORM = 15.1581923699 
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18) 


19) 
29) 


21) 


22) 


23) 


24) 


Comments 


Run program again. 


Convert from rectangular 
to cylindrical. 


Change output format to 
match conversion. 


Specify rectangular input. 
Specify cylindrical output. 


Choose conversion now that 
I/O format is correct. 


Input vector. 
Display result. 


End program. 


EXAMPLES 


<any key> 
[Y] 


[C]} 


[R] 
[C] 
[Cc] 
-4 [ENDLINE] 


7 [ENDLINE] 
@ [ENDLINE] 


<any key> 
<any key> 


<any key> 
[N] 
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Display 


RUN AGAIN? (Y/N) 
A,S,X,1,D,N,M,U,C,F,Q? 
I/O INCORRECT 

INPUT R/C? 

DISPLAY R/C? 
A,S,X,1,N,M,U,C,F,Q? 
x= 


y= 
Z= 


r= 8.0622577483 
theta= 119.744881297 
z= @ 


RUN AGAIN? (Y/N) 
> 


0819 
8828 
8036 
0849 
8858 


LISTING 


! VECTOR OPERATIONS 
! Revision 1.8 3/28/84 
t 


DEGREES @ DELAY @,@ @ CFLAG 1,2,3,4 

K@S="ASXIDNMUCFQ" @ I@S="CR" @ I1S="YN" 

'IO': 

DISP "INPUT R/C ? "s 

'I1l': K1S=UPRCS (KEYS) @ IF NOT POS(I@$,K1$) THEN 'I1' ELSE DISP 
IF K1S$="C" THEN SFLAG 1 ELSE CFLAG 1 

IF FLAG(1) THEN X1S="r="_@ Y1S="theta=" ELSE X1S="x=" @ Y1lS="y=" 
DISP "DISPLAY R/C ? "; 

'Ol': K1S=UPRCS (KEYS) @ IF NOT POS(I@S$,K1$) THEN 'O1' ELSE DISP 
IF K1$="C" THEN SFLAG 2 ELSE CFLAG 2 

IF FLAG(2) THEN X2S="r=""@ y2S="theta=" ELSE X2S="x=" @ Y2S="y=" 
"OP's: DISP "A,S,X,1,D,N,M,U,C,F,Q" 

'WAIT': KIS=UPRCS(KEYS) @ K@=POS(K@$,K1S) @ IF NOT K@ THEN 'WAIT' 
IF K@=11 THEN 'Q' 

IF K@=10 THEN CFLAG 1,2,3 @ GOTO 'IO' . 
IF KO#9 THEN 'Al' ae 
IF NOT FLAG(1) EXOR FLAG(2) THEN DISP 'I/O INCORRECT' @ WAIT 2” Al‘ 
GoTo 'IOo' 

‘Al's: IF K@<6 THEN 'I2VA' 

‘Ilv': IF FLAG(3) THEN 'SUBCALL' 

DISP X1$; @ INPUT "";X1 

IF FLAG(1) AND X1<@ THEN DISP "INVALID ENTRY" @ GOTO ‘Ilv' 

DISP Y1$; @ INPUT "";Y1l 

INPUT “z="3Z1 

GOTO 'SUBCALL' 

'I2VA': IF FLAG(3) THEN 'I2VB' 

DISP “VECTOR 1:"&X1$; @ INPUT "";X1 

IF FLAG(1) AND X1<@ THEN DISP “INVALID ENTRY" @ GOTO 'I2VA' 

DISP "VECTOR 1:"&Y1S; @ INPUT "";Y1 

INPUT "VECTOR 1:z="321 

"T2VB's: 

DISP "VECTOR 2:"&X1S; @ INPUT ""*x2 

IF FLAG(1) AND X2<@ THEN DISP "INVALID ENTRY" @ GOTO 'I2VB' 

DISP "VECTOR 2:"&Y1S; @ INPUT “""sy¥2 

INPUT "VECTOR 2:z2=";22 

IF FLAG(1) THEN CALL C2R(X2,Y2) 


"SUBCALL': 

IF FLAG(1) THEN CALL C2R(X1,Y1) 
GOSUB K1$ 

‘CYCLE’: 

DISP "RUN AGAIN? (Y/N) "; 


'A2's: K1LS=UPRCS(KEYS) @ ON POS(I1S,K1S)+1 GOTO 'A2','USEIT','Q' 
‘USEIT': DISP @ IF IP((K@-1)/3)=l OR K@=9 THEN CFLAG 3 @ GOTO ‘OP' 
DISP “USE RESULT? (Y/N)"; 

'A3': KIS=UPRCS(KEYS) @ IF NOT POS(I1$,K1$) THEN 'A3' ELSE DISP 
IF POS("N",K1S) THEN CFLAG 3 @ GOTO 'OP' 

IF NOT POS("Y",K1S) THEN 'USEIT! 

IF FLAG(2) THEN CALL €2R(X@,YQ@) 

X1=xX@ @ Yl=Y@ @ Z1=Z20 @ SFLAG 3 @ GOTO ‘OP' 

ASS 

XO=X1+X2 @ YO=Y1+Y2 @ ZG=Z1+Z2 

GOSUB ‘'OV' 
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1858 
1969 
1879 
1880 
1699 
1168 
1116 


LISTING 


RETURN 
"Ss 

X@=X1-X2 @ Y@=Y1-Y2 @ 720=71-22 

GOSUB ‘'OvV' 

RETURN 

[DY 

DO=X1*X2+Y1*Y2+Z1*2Z2 

DISP “DOT PRODUCT =";D8 @ GOSUB ‘WwW! 
RETURN 

"x6 3 

XO=Y1*Z2-Y2*Z1 @ Y@=Z21*X2-Z2*X1 @ ZO=X1*Y2-xX2*Y1 
GOSUB '‘Ov' 

RETURN 

‘M's: INPUT "SCALAR ="s5S 

X@=X1*S @ YO=Y1*S @ Z2ZG=Z1*S 

GOSUB 'ov' 

RETURN 

hose 

SFLAG 4 @ GOSUB 'N' @ CFLAG 4 

M3=SQR ( (X2-X1) *2+(Y2-yY1) *2+ (22-21) *2) 
A@= (M1*M14+M2*M2-M3¥*M3) /(2*M1*M2) 
C9=.G9GAGHBAGHBG2 

IF FP (FP (ABS (AS) )+C9)<. 000090000045 THEN AG=IP (AG+SGN (AQ) *C9) 
A1=ACOS (AG) 

DISP “INCLUDED ANGLE =";Al @ GOSUB 'wW' 
RETURN 

'n's 

M1=SQR(X1*X1+Y1*Y1+21*Z1) 
M2=SQR(X2*X2+Y2*Y2+Z2*Z2) 

IF NOT FLAG(4) THEN DISP "NORM ="s:M1 @ GOSUB 'w' 
RETURN 

‘u's 

SFLAG 4 @ GOSUB 'N' @ CFLAG 4 

XO=X1/M1 @ Y@=Y1/M1 @ 20=21/M1 


GOSUB 'oOvV' 

RETURN 

‘W's: IF KEYS="" THEN 'W' ELSE RETURN 
"ers 

X@=X1 @ Y@=Y1 @ ZG=Z1 

‘Ov': 


IF FLAG(2) THEN CALL R2C(X@,Y@) 
DISP X2$;X@ @ GOSUB ‘w' 
DISP Y2S;Y@ @ GOSUB '‘'W' 
DISP "z=";Z8 @ GOSUB 'W' 
RETURN 
'Q': CFLAG 1,2,3,4 @ DISP @ PUT "#43" @ END 
SUB R2C(R,T) 
I=R @ J=T 
R=SQR(I*I+J*J) 
T=ANGLE (I,J) 
SUB END 
SUB C2R(I,J) 
R=I @ T=J 
I=R*COS (T) 
J=R*SIN(T) 
SUB END 
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Numerical Integration 


« This program will perform numerical integration whether a function 
is known explicitly or at a finite number of equally spaced points. 








For the explicit case, a 16 point Gaussian quadrature is provided 
for a finite interval. Also, the integral may be found for the 
explicit case using Simpson's or the Trapezoidal methods. The 
quadrature is more accurate, however. 


Given a finite set of points, Simpson's rule or the Trapezoidal 
rule may be used to find the integral. 


Equations: 
Gaussian Quadrature: 
N 
(b-a)/ 2) YY wk<j>f( (atb+x<4> (b-a) ) /2) 
jal 


where w<j>,x<j> are the weights and nodes, respectively. The 
weights and nodes are for an integral from -1 to 1. 


Simpson's Rule: 


h[f£(x@) + 4£(x1) + 2£(x2) + .. + 2£(x<n-2>) + 4£(x<n-1>) + £(x<n>)]/3 


( where h is the (equally spaced) interval distance. 
‘n' must be even for the method to work. 


Trapezoidal Rule: 


n=l 

h/2[£(x9) +2 DO (£(x<j>) + £(x<nd)] 
j=l 

where h = (b-a)/n for interval (a,b). 
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User Instructions 


Comments Input Display 
J --- a nn rere T-------------- T2222 nn ner I 
1) Run the program. Gauss,Simp,Trap,Quit 
2) Choose desired method. 
a. Gauss G f(x) = 
b. Simpson Ss Explicit (Y/N) 
c. Trapezoidal T Explicit (Y/N) 
d. Quit Q <blinking cursor> 
Go to desired option to 
continue. 
kkkkkk Gauss kkkkke 
Gl) Key in the function with 
variable 'x',. <function> 
[ENDLINE] Lwr,Upr bounds= 
G2) Key in the lower and upper 
bound separated by a comma. <L,U> calculating.. 
Result= <val.> 
G3) To continue, press any key. <any key> Gauss,Simp,Trap,Quit 
Go to step 1 to continue. 
kkkkKKK Simpson kkekekekee 
Sl) If the function is known: Y # of Partitons=<val.> 
Go to the explicit case (S5). 
If the function is unknown: N Nmbr of pnts=<val.> 
$2) (Implicit case) 
Key in the number of points. 
For Simpson's rule, this must 
be an odd number. <number> Interval length=<val.> 
$3) Key in the interval length 
for the partitions. This 
routine requires that parti- 
tions be the same length. <length> E (0) =<val.> 
You are now in the matrix 
editor. Continue at that 
section (El) and return. calculating.. 
Result=<val.> 
S4) To continue, press any key. <any key> Gauss,Simp,Trap,Quit 
Go to step 1 to continue. 
$5) (Explicit Case) 
Key in the number of 
Partitions desired. This must 
be an even value. <number> £(x)= 
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User Instructions 


Comments Input Display 
Tm a nr ee [-------------- Ta en ne ne ee I 
S6) Key in the function with 
variable 'x', <function> 
[ENDLINE] calculating... 
Result= <val.> 
S7) Key in the lower and 
upper bound separated by 
a comma. <lb,ub> calculating... 
Result= <val.> 
S8) To continue, press any key. <any key> Gauss,Simp,Trap,Quit 


Go to step 1 to continue. 


kak kkk Trapezoidal kRekKEKE 


Tl) 


T2) 


T3) 


T4) 


T5) 


T6) 


T7) 


T8) 


If the function is known: 
Go to the explicit case (T5) 
If the function is unknown: 


(Implicit case) 
Key in the number of points. 


Key in the interval length 
for the partitions. 


You are now in the matrix 
editor. Continue at that 
section (El) and return. 


To continue, press any key. 
Go to step 1 to continue. 


(Explicit Case) 
Key in the number of 
partitions desired. 


Key in the function with 
variable 'x'. 


Key in the lower and 
upper bound separated by 
a comma. 


To continue, press any key. 
Go to step 1 to continue. 


Y # of Partitons=<val.> 

N Nmbr of pnts=<val.> 
<number> Interval length=<val.> 
<length> E(@)=<val.> 


calculating.. 
Result=<val.> 


<any key> Gauss,Simp,Trap,Quit 


<number> £ (x)= 


<function> 
{ENDLINE] calculating.. 
Result= <val.> 


<1lb,ub> calculating... 
Result= <val.> 


<any key> Gauss,Simp,Trap,Quit 
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User Instructions 


Comments Input Display 
I------------~------+-----------. I-------------- [--~---------~-------~- I 


REKEKE Matrix Editor *****x*x 


The points E(@) through E(n) 

May be entered. The value 
displayed is the current value 
of the point. You may change the 
value, move to a previous value, 
move to the next value, or 
specify a point to move to. 
Pressing [Q] exits the matrix 
editor, and the program finds 


the value of the integral. E (I) =<value> 
El) To change value: [ENDLINE] 

<new val.> 

{ENDLINE] E(I+1)=<value> 


To move to next 

value (this is the 

arrow key, not the 

greater than char.): {>] E(I+1) =<value> 


To move to previous 
value (this is the 
arrow key, not the 


less than char.): {<] E(I-1) =<value> 

TO move to a specific 

point: [SPC] Element= 

Key in the desired 

element index: <index> E(<index>) =<val.> 
To quit the routine: {Q] calculating... 
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EXAMPLE PROBLEMS 
Problem 1. 
Given the approximations below for £(x), compute the approximations 


to the integral from the bounds @-2 by the trapezoidal rule and 
by Simpson's rule. The interval length is 9.25. 


1 G 1 2 3 4 5 6 7 8 

f£ (x) 2 2.8 oe 562 7 9.2 12.1 15.6 20 
Comments Input Display 

I-----------------~------------- Pase445e-o esse J----------------+-+-~ I 

1) Run the program. Gauss,Simp,Trap,Quit 


2) Choose trapezoidal 
method: T Explicit (Y/N) 


3) Since the function is 
unknown, choose implicit 


case: N Nmbr of pnts=@ 
4) Key in the number of points. Qfexoiwe} Interval length=@ 
5) Key in the interval length 
for the partitions. o25feworwa] E(B) =G 
6) Enter the function values: [ENDLINE] value= 
2 
[ENDLINE] E(1)=@ 
[ENDLINE] value= 
2.9 
[ENDLINE] E(2)=@ 


The previous value was 
supposed to have been 2.8. 


Edit it. [<] E(1)=2.9 
[ENDLINE] value = 
2.8 
[ENDLINE] E(2) =@ 


Continue entering the rest of 
of the values in the same manner. 


To quit: Q calculating.. 
Result=16.6759@ 
7) To continue, press any key. <any key> Gauss,Simp,Trap,Quit 
8) Now use Simpson's rule. Ss Explicit (Y/N) 


9) Since the function is 
unknown, choose implicit 
case: N Nmbr of pnts=9 
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User Instructions 


Comments Input Display 
[----~-------~------------------- I-------------- I---~~----~------------ I 
18) Since this is the correct 

number, just press: [ENDLINE] Interval length=.258@ 
11) This value is also okay. [ENDLINE] E(@)=2 


12) Since the program leaves 
the values the same, simply 


press [Q] to finish. Q calculating... 
Result=16.5833 
13) To continue, press any key. <any key> Gauss,Simp,Trap,Quit 
Problem 2. 


Find the value of the integral of f(x) from @ to 2*pi where 
£(x) = 1/(1l-cos(x)+ .25). 


Comments Input Display 
I----------------~-------.~------ I-------------- [--------~--~-------~---- I 
1) Run the program. Gauss,Simp,Trap,Quit 
2) Choose Gauss's 

method: G £(x)= 
3) Key in the function with 
variable 'x', 1/ (l-cos (x)+.25) 
[ENDLINE] Lwr,Upr bounds= 


4) Key in the lower and 
upper bound separated by 


a comma. 0@,2*pi 
[ENDLINE] calculating.. 
result= 8.3776 
5) To continue, press any key. <any key> Gauss,Simp,Trap,Quit 
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LISTING 
Note that some of the comments are not preceded by line numbers. 


9919 DEF FNKS(D@S,K@S) ! Find which key was pressed function 
@028 DISP D@S 

0030 'KEY': KS=KEYS 

0948 IF NOT POS(K@S,KS) THEN "KEY" 

0958 FNKS=KS 

00680 END DEF 


RakKKKKKKKRKRKKKK INITIALIZATION *¥R#REKKKKKER 


6876 OPTION BASE @ 

8888 DESTROY A,H,P,TS$,KS,K1$ 
8698 DIM K$[4],K1$[1],TS$[1] 
0190 INTEGER N,P 


kkk*k BEGIN **** 


01190 ‘START': K1S=FNKS ("Gauss,Simp,Trap,Quit","GSTQ") 
9128 GOTO K1S 


*x*k* SIMPSON, TRAPEZOIDAL START **** 


9130 'S': 'T':s TS=FNKS("Explicit?(Y/N)","YN") 
9149 IF TS="Y" THEN "PRT" 


**e*k* For implicit case, get function values 


8150 ‘PL': INPUT “Nmbr of pnts= ",STRS(N);N 

0160 IF N<=@ THEN DISP "Must be positive" @ GOTO 'PL' 

0179 IF K1S="S" AND MOD(N-1,2) THEN DISP “MUST BE ODD NUMBER" @ GOTO 
"w PL " 

186 INPUT "Interval Length= ",STRS(H) ;H 

0198 DIM A(N-1) 

0200 CALL EDIT1(A,N-1) 

6218 GOTO "CALC" 


xe*k* For explicit case, get number of partitions **** 

9220 '"PRT': INPUT "# of Partitions= ",STRS(P);P 

@230 IF P<=@ THEN DISP “Must be positive" @ GOTO 'PRT' 
8248 IF K1$="S" AND MOD(P,2) THEN DISP “MUST BE EVEN, NONZERO" @ GOTO 
Ww PRr" 

*#*k*k* For explicit case, get function **** 

@250 'G': INPUT "f(x)= ",FS;FS 

0260 IF K1$="G" THEN TS="N" 

0270 IF FS="" THEN "G" 

****k Get the bounds for the integral **** 

$8289 'B': INPUT "Lwr,Upr bounds= ",STRS(L) &","&STRS(U) ;L,U 
6298 IF L=U THEN DISP “THE INTEGRAL VALUE IS @" @ GOTO "B" 
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LISTING 


$300 'CALC': O=FLAG(-16,0) @ T=FLAG(-19,1) ! Temporarily set radians 
6318 DISP “calculating.." 


8328 GOSUB 


K1S&TS$ 


0330 O=FLAG(-16,0) @ T=FLAG(-198,T) ! Re-establish angular mode 
EKS ("2F946",4) @ DELAY INF,INF ! Save wait period 


8348 D1S=PE 
8358 DISP " 
8368 POKE " 
8378 GOTO " 


kKkkKKK SUB 


Result=";R 
2F946",D15 
START" 


PROGRAM CALL 


Display result 
Re-establish wait period 


SECTION ****** 


8388 'GN': CALL GAUSS(FS,L,U,R) @ RETURN 

0398 ‘SN': CALL SIMPSON (A,N-1,H,R) @ RETURN 

9490 'TN': CALL TRAP(A,N-1,H,R) @ RETURN 

8410 ‘SyY's: CALL SIMPSONE(FS,L,U,(P),R) @ RETURN 

420 ‘Ty's: CALL TRAPE(FS,L,U,(P),R) @ RETURN 

9430 'Q': PUT "#38" @ END ! Restore cursor and end program 
044¢ ! 


**x*e* Trapezoidal rule, explicit case **** 


8458 SUB TRAPE(FS,L,U,P,R) 
0468 H=(U-L)/P 


8476 Q=9 
9488 X=L @ 
6498 X=U @ 


9580 FOR I= 


8518 X=L+H* 


9530 NEXT I 
9548 R=(R+2 


R=VAL (FS) 
R=R+VAL (FS) 
1 TO P-1l 


I ! what is A? 
@528 Q=Q+VAL (FS) 


*Q) *H/2 


6556 END SUB 


0568 ! 


**ek* Trapezoidal rule, 


implicit case **** 


9578 SUB TRAP(A(),N,H,R) 


0580 R=A(@) 
9598 Q=98 


9600 FOR I= 


8618 Q=Q+A( 
66208 NEXT I 
9630 R=(R+2 


+A (NN): 


1 TO N-1 
T) 


*Q) *H/2 


9648 END SUB 


6658 !} 


**** Simpson's rule, 


implicit case **** 


8668 SUB SIMPSON(A(),N,H,R) 


0670 R=A(@) 


+A (N) 


968@ FOR I=1 TO N-3 STEP 2 
8690 R=R+4*A (I) +2*A(I+1) 


$700 NEXT I 
8719 R=(R+4 


*A(N-1)) *H/3 
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( 


8720 
873G 


kkkk 


8749 
8758 
87698 
6778 
8780 
8798 
8888 
88190 
8828 
6830 
8840 
8859 


kik 


6860 
8870 
8889 
8898 
8988 
89198 
8928 
8936 
G94 
8950 
8960 
89708 
8986 
8990 
1988 
1616 
1920 


kRkKK 


103¢ 
1848 
1859 
1968 
1878 
1989 
1699 


1168 


1119 
112g 
1139 
1149 


LISTING 


END SUB 
{ 


Simpson's rule, explicit case **** 


SUB SIMPSONE(FS,L,U,P,R) 
H=(U-L) /P 

X=L @ R=VAL(FS) 

X=U @ R=R+VAL (FS) 

FOR I=l TO P-3 STEP 2 
X=L+H*I @ R=R+4*VAL (FS) 
X=X+H @ R=R+2*VAL (FS) 
NEXT I 

X=U-H 

R= (R+4*VAL (FS) ) *H/3 
END SUB 

! 


16 point Gaussian method **** 


SUB GAUSS (FS,L,U,R) 

DIM N(1,7) 

DATA .989466934992,.944575923073,.865631202388,.755404498355 
DATA .617876244403,.458016777657,.281603550779,.950125098376E-1 
DATA .271524594118E-1,.622535239386E-1,.951585116825E-1,.124628971256 
DATA .149595988817,.169156519395,.182603415045,.189450610455 
READ N(,) 

R= 

C=(U-L) /2 

FOR I=@ TO 7 

X=N(@,1I) *C+(L+U) /2 

R=N(1,1) *VAL(FS)+R 

X=-N (0,1) *C+(L+U) /2 

R=N(1,1) *VAL(FS$)+R 

NEXT I 

R=R*C 

END SUB 


Matrix Editor **** 


SUB EDIT1(A() ,N) 

DEF FNK$(D@$,K@$) ! Find which key was pressed function 
DISP D@S 

"KEY': KS=KEYS 

IF NOT POS(K@$,KS$) THEN "KEY" 

FNKS$=K$ 

END DEF 


DEF FNFS(Y) ! Fix display function. Temporarily alters, 
! then returns to original setting. 

DS=PEEKS ("2F6DC",2) @ STD 

ENFS=STRS (Y) 

POKE "2F6DC",DS 

END DEF 


nada 19 


1158 
1166 
1176 
1180 
1198 
1266 
1218 
1226 
123 
1246 
1256 
1260 
1278 


1286 
1298 
1380 


1316 
1320 


LISTING 


I= 

"LOOP': K1S=FNKS ("E("&FNFS (I) &") = "&STRS(A(I)),"Q #38#47#48") 
IF K1S="_" THEN 'M'! 

IF K1S="Q" THEN 'Q! 

IF LEN(K1$) #3 THEN 'LOOP' 

K1$[1,1]="k" 

GOSUB K1S$ 

GOTO 'LOOP' 

‘K47': I=MOD(I-1,N+t1) ! Move to previous position 
RETURN 

"K38': INPUT "“Value= ";A(I) 
*K48': I=MOD(I+1,N+1) 
RETURN 


Change value 
Move forward 


o— oe 


‘M's: INPUT "Element= ";I ! Get desired position 
IF I<@ OR I>N THEN DISP “EXCEEDS BOUNDS" @ WAIT 1 @ GOTO 'M' 
GOTO "LOOP" 


‘o's PuT "#38" ! Quit routine 
END SUB 
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Solution to F(x) = 0 on an Interval 








This program provides two methods to find a real root of the 
equation f(x) = @. They are Newton's Method and the Pegasus Method. 
In addition, the program allows the user to find the value of the 
function for an input x. 


Input for the Newton's method consists of the function to be solved, 
one initial guess, and as an option, the derivative of the function. 
If the derivative is not input, a numerical approximation is used. 


Input for the Pegasus method consists of the function to be solved 

and two initial guesses that must bound the root. This implies that 
£(x@)*£(x1)<@. The routine to calculate function values may be used 
to establish a legal interval. 


When a root is found, the output will consist of the x value 
displayed to the setting of the computer and the value f(x), where x 
is the displayed root. It is possible that f(x) will not be exactly 
®@. However, it will generally be within an acceptable range around 
zero based on the number of significant digits in the input. If the 
desired accuracy is not obtained, it may be possible to decrease the 
value used to check for acceptance (the variable E in both sub 
programs). In some instances, the function may have to be modified. 


Newton's method converges to a root quickly in cases where it can 
find one. Its ability to locate a root depends on the function and 
the initial guess. It is not guaranteed to find a root. If the 
derivative is @ or 5@ iterations are peformed, the routine exits, 
displaying an appropriate message to the user. The number of 
iterations may be changed by altering the program. 
The Pegasus Method is a modified regula falsi method with an 
estimated order of convergence superior to a secant method. For any 
legal interval, the method is guaranteed to converge. 
The equations: 
Newton's Method: 

x<nt+l> = x<nd>-f(x<n>) /£' (x<n>) 

The exit criteria is abs(x{n+l]-x[{n]) <= epsilon 

where epsilon is a small value. 


£'(x) Approximation: 


When the derivative is not given, the program uses 
the following approximation: 


E'(x) <== (£(x+1I/2) - £(x-1/2))/1 


where I = .80@1(x) or .@G0001 if x= @. 
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Pegasus Method: 
The Regula Falsi method used is: 
x<n+l> = x<n> - £(x<n>) [ (x<n>-x<n-1>) /(£ (x<n>) -f (x<n-1>))] 
The approximations for the next iteration are chosen by: 


if £(x<n+l>)£(x<n>) < @ then 
(x<n-1>,f£(x<n-1>)) <== (x<n>, £(x<n>)) 
if f(x<n+1>)£(x<n>) > @ then 
(x<n-1>,£(x<n-1>)) <== (x<n-1>, f£(x<n-1>) /(f (x<n>) +f (x<n+1>))) 


References: 


Dowell, M. & Jarrett, P.; "The Pegasus Method for Computing the Root 
of an Equation", BIT 12 (1972) pp. 503-568 


Atkinson, Kendall E., "An Introduction to Numerical Analysis", Wiley 
and Sons 


Carnahan, B., Luther, H.A., and Wilkes, J.0., “Applied Numerical 
Methods", Wiley and Sons, Inc. 
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USER INSTRUCTIONS 


Comments Input Display 
[---~--------------------------- I-------------- [--------------------- I 
1) Run program. £(x)= 
2) Key in the function using <function> Root, F(x) ,Chngf,Quit 


the character 'x' as the 
variable. 


3) Press the capital letter 
of the desired operation. 


'R' begins the solve R Pegasus, Newton 
routine. 

‘F' finds function F X= 
values for given x's. 

'C' allows the function Cc £ (x) =<function> 
to be changed. 

'Q' quits the program. Q Done 


Go to the appropriate 
heading to continue. 


REEKEX ROOT OPTION #8 RRR Pegasus, Newton 
Rl) Choose desired method by 


pressing the appropriate 
capital letter. 


For Pegasus, press P. P lower bound: 
Go to step R2. 
For Newton, press N. N derivative= 


Go to step R8. 


R2) PEGASUS METHOD lower bound: 
R3a) Key in a lower bound. <value> 

{ENDLINE] upper bound: 
R3b) To exit, press: [ENDLINE] Root,F(x) ,Chngf,Quit 
R4a) Key in upper bound. <value> 

[ENDLINE] calculating.. 
R4b) To exit, press: [ENDLINE] Root,F(x) ,Chngf,Quit 
R5a) The answer will be x= <result> 

displayed. 


R5b) If the interval does 
not bound a root a 
message will be dis- 
played and you will 
return to step R4. 


R6) To see the function 


value at the root, 
press any key. <any key> f (x) =<value> 
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USER INSTRUCTIONS 


Comments Input Display 
[----------~-------------------- [-------------- I--------------------- I 
R7) To exit, press any key. <any key> Root,F (x) ,Chngf,Quit 
R8) NEWTON METHOD derivative = 
R9a) If you want to use the 

derivative, key it in. <derivative> 
[ENDLINE] initial guess= 
R9b) If you do not wish to 
enter the derivative: [ENDLINE] initial guess= 
R108) Key in initial guess. <guess> 
[ENDLINE] disp convergence? 
Rll) If you wish to watch 
the convergence, press Y or calculating.. 
‘y', else press any <another key> 
other key. 
R12) If a root is found, the x= <result> 


result will be found. 
If a root is not found, 
an appropriate error 
message will be given. 


R13) To see the function 
value at the root, 


press any key. <any key> f (x) =<value> 
R14) To exit, press any key. <any key> Root,F(x) ,Chngf,Quit 
xekkkEk FIND FUNCTION VALUES ****** X= 
Fl) Key in the value of x <value> 
at which you wish the [ENDLINE] f(x) = <result> 
function evaluated at. 
F2) Press any key to <any key> X= 
continue. 


F3a) Continue at step F2 for 
more values. 


F3b) To exit, press: [ENDLINE] Root, F(x) ,Chngf,Quit 
*#eeeEX CHANGE FUNCTION ****** £(x)=<crnt function> 
Cl) Key in desired function. <function> 

{ENDLINE] Root,F(x),Chngf,Quit 
*xxkkkE QUIT PROGRAM ****** Done 
Ql) To get the prompt back: <any key> <blinking prompt> 


press any key. 
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EXAMPLE 


This example demonstrates the various options of the program. It 
must be started at the beginning and followed to completion for the 
given keystrokes to work as listed. The display should be FIX 1l. 


Find a root for the function f(x) = In(x) + 3x - 10.8074. Use the 
Pegasus method, then Newton's method. Then find a root for the 
function £(x) = 3x*6 - 22x75 +11. 


Comments 


Display 


[----=----------------.--.------ feesbees eee cee I-----------------.--- I 


1) Run program. 


2) Key in function. Note that 
it is entered in a BASIC 
format. 


3) Call the solve section. 

4) Choose Pegasus method. 

5) Guess at a bound: 

6) Guess at an upper bound. 
7) The interval did not bound 


the root. Use F(x) to find 
an interval. 


8) Find £(5). 


9) Find £(190) 


186) Try £(1) 


ll) Since the function is 
continuous on the 
interval [1,5], and the 
values of the function 
at these points are of 
opposite sign, this 
interval bounds a root. 
Exit and move to the 
solve section. 


12) Key in lower bound 


In (x) +3*x-16.8074 


[ENDLINE] 


R 
P 
5 [ENDLINE] 


18 [ENDLINE] 


[ENDLINE] 

F 
5 [ENDLINE] 
[ENDLINE] 


1@ [{ENDLINE] 
[ENDLINE] 


1 [ENDLINE] 
[ENDLINE] 


[ENDLINE] 
R 
P 


1 [ENDLINE] 


page 25 


Root,F(x),Chngf,Quit 


Pegasus ,Newton 
lower bound: 
upper bound: 


intrvl must bound root 
lower bound: 


Root, F(x) ,Chngf,Quit 
x= 


£(x) = 5.8020379124 
x= 


21.4951856893 


x= 


£(x)= -7.8074 
x= 


Root,F (x) ,Chngf,Quit 
Pegasus,Newton 
lower bound: 


upper bound: 


[-------- ~~ + + + - - - - - - -- - = 


14) 


15) 


16) 


17) 


18) 
19) 


20) 


21) 


22) 


Comments 
Key in upper bound. 
See how close £(3.21..) 
is to @. 
Continue 
Solve using Newton's 
Method. 
Let routine approximate 
the derivative. 
Key in an initial guess. 


Let's not display it. 


We must now solve the 
second function. 


Use the derivative 
this time. 


The intermediate 
results are displayed. 


Quit the program. 


EXAMPLE 


5 [ENDLINE] 


[ENDLINE] 


{ENDLINE] 


{ENDLINE] 

1 [ENDLINE] 
N 

{ENDLINE] 

[ENDLINE] 


Cc 
3*x*6-22*x°54+11 

R 

N 


18*x*5-110*x%4 
{ENDLINE] 
5 [ENDLINE] 

Y 


[ENDLINE] 
[ENDLINE] 


Q 
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Display 


calculating.. 
X= 3.21336987018 


f(x) = @ 


Root,F (x) ,Chgf,Quit 


Pegasus, Newton 
derivative= 


initial guess= 
disp convergence? 


calculating.. 

X= 3.21336087017 
f(x)= @ 

Root,F(x) ,Chngf ,Quit 


f (x) =LN(X)+3*X-10.8074 
Root,F (x) ,Chngf ,Quit 
Pegasus ,Newton 
derivative= 


initial guess= 
disp convergence? 
calculating.. 
2.25088 
2.4794542177 
1.93158885068 


- 89346784412 

x= .893467844031 
£(x)= @ 

Root,F(x) ,Chngf ,Quit 


DONE 


EME n 


LISTING 
Note that some of the comments are not preceded by a line number. 


@018 ! SOLUTION TO F(X)=@ 

6926 ! REV 1.8 -- 1/25/84 

0030 FS="" 

0040 ‘C's: INPUT "f£(x)=",FS;FS ! get function 

6858 'D': DISP “Root,F(x) ,Chngf,Quit" ! main prompt 

@@6@ AS=KEYS @ IF AS="" THEN 69 

0070 IF POS("RFQC",UPRCS(AS[1,1])) THEN GOTO UPRCS(AS[1,1]) ELSE GOTO 'D' 


kkkKKKK Pind the Root kkkRKKK 


0080 'R': DISP "Pegasus,Newton" ! get desired method 

0890 MS=KEYS @ IF MS="" THEN 4@. 

0190 IF POS("PN",UPRCS(M$[{1,1])) THEN GOSUB M$[1,1] ELSE GOTO "Rr 
0110 IF R<>INF THEN GOTO ‘'DISPR' 

120 DISP "NO ROOT FOUND" 

81308 GOTO ‘'RTN' 

0149 'DISPR': X=R ! display result 

@15@ DISP "x= "3x 

0160 AS=KEYS @ IF AS="" THEN 164 

0176 DISP "£(x)= ";VAL(FS) ! display function value 

@18@ 'RTN': AS=KEYS @ IF AS="" THEN GOTO 'RTN' ELSE GOTO 'D' 


ekkkkKK Find values of £(x) given x ******* 


6190 'F':.ON ERROR GOTO 'XIT' 

G20 INPUT "X= ";X 

$210 Y=VAL(FS) 

G@228 DISP "f(x)= ";Y 

9230 AS=KEYS @ IF AS="" THEN 236 ELSG@ “F" 
8248 GOTO 'F' 

@258 'XIT': OFF ERROR @ GOTO 'D' 


eeeKEEKK Set up to call Pegasus method ******* 


@26@ 'P': ON ERROR GOTO 'XIT' 

8276 INPUT "lower bound:";L 

0280 INPUT “upper bound: ";U 

9290 X=L ! see if interval contains root 

0308 Y1=VAL (FS) 

0318 X=U 

6320 Y2=VAL(ES) 

0330 IF Y1*Y2>@ THEN DISP “intrvl must bound root" @ GOTO 'P' 
0340 CALL PEG(FS,L,U,R) 

9358 RETURN 


kekKKEK Set up to call Newton method ******* 


9368 'N': DESTROY NS 

@370 ON ERROR GOSUB 'NEWTERR' ! catches no derivative option 
6388 INPUT “derivative= ";DS$ 

0390 INPUT “initial guess= ","1";xX@ 

8480 DISP "disp convergence?" 
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LISTING 


@41@ AS=KEYS @ IF AS="" THEN 4340 

0420 IF UPRCS(AS$)="Y" THEN D=1 ELSE D=@ 
9438 CALL NEWT(FS,DS,X@,D,R) 

9449 RETURN 

0459 ‘NEWTERR': OFF ERROR 

8460 ON ERROR GOTO 'XIT' 


kkeekkek Quit program kkkekekeek 
0470 'Q': DISP “™ DONE" @ END 


kkkkKKK Pegasus subprogram ******* 

The inputs are the function (FS), lower limit (X@), and upper limit 
(Xl). The result is returned through variable 'R'. The function 
string must be a legal BASIC expression in variable 'x'; the limits 
must bound a root. 


@480 SUB PEG(FS,X%,X1,R) 

@490 DISP “calculating.." 

0500 E=.@000000001 ! error tolerance 
0510 C=2*X1-xd 

0520 X=X@ ! init y@,yl 

0538 Y@=VAL (ES) 


8549 X=X1 
8558 Y1=VAL (FS) 
8560 ! 


6570 'LP’: X2=X1-Y1* ((X1-xX@) /(Y1-Y@) ) 
0588 IF ABS(X2-C)<=E THEN R=X2 ELSE GOTO ‘CT! 
8598 GOTO 'ND' 

606 'CT': X=X2 

8618 Y2=VAL (FS) 

8628 Cl=Y2*Y1 

6630 IF Cl<® THEN X@=X1 @ Y9=Y1 

6649 IF Cl>@ THEN YO=YO*Y1/(Y1+Y2) 
8658 X1=X2 

@668 Y1=Y2 

8678 C=X2 

688 GOTO ‘LP’ 

8698 'ND': SUB END 


week Newton subprogram ******* 

The inputs are the function (FS), optional derivative (D$), initial 
guess, and the display-convergence boolean. The result is returned 
through 'R'. The functions must be in variable 'x'. If no root is 
found, 'R' is set to the value ‘inf'. 


9790 SUB NEWT (FS,DS,X@,D,R) 

8716 DISP “calculating.." 

0720 B=.0000000001 ! error tolerance 

0730 L=8 ! INIT LOOP COUNTER 

0740 'LP': IF DS="" THEN GOSUB 'AD' ELSE X=X@ @ Y1=VAL(DS) 
0750 IF Yl=@ THEN R=INF @ DISP “DERIVATIVE=9" @ GOTO 'ND' 
0768 X=x@ 

0778 Y=VAL(FS) 
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8788 
8798 
8808 
0818 
8828 
8830 
8848 
6850 
8868 
6878 
8880 
8890 
8988 
G91@ 


X1=x-yY/Y1 
IF ABS (X1-X)<=E THEN R=Xl1 @ GOTO 'ND' 


X8=X1 
L=L+1l 


IF D THEN DISP X@ 


LISTING 


IF L=5@ THEN DISP "5@ ITERATIONS" @ R=INF @ GOTO 'ND' 
GOTO 'LP' 


‘AD': 


IF X@=0 THEN I=.000001 ELSE I=.0901*xX@ 


X=XO+1/2 


Y1=VAL 


(FS) 


X=XG-1/2 


Yl=(Y1-VAL(FS$))/I 


RETURN 
'ND': 


SUB END 
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deriv. 


approx. 


Matrix Operations 


eee SSS SSS. ST ewww 


This program allows the user to caculate the determinant of a real 
valued matrix and find the inverse or solve a system of equations 
for real or complex valued systems. 


The method used is Gaussian elimination with partial pivoting. The 
matrix is decomposed into an LU form and the pivoting strategy is 
kept track of in a separate matrix. For a more in-depth discussion, 
please refer to the references. 


The determinant is calculated from the decomposed matrix by 
multiplying the values in the main diagonal. It is found only for a 
real valued matrix. 

The inverse is found by solving the system Ax(i) = I(i), where T(i) 
is the ith column of the identity matrix and x(i) is the ith column 
of A inverse. This is performed N times as i ranges from 1 to N. 

If the inverse or solution to a system of equations is attempted that 
involves a singular matrix, the message "MATRIX IS SINGULAR" will be 
displayed. 

Remarks: 


The program will use the particular number display you specify 
before running the program. 


If you pause the program during the listing of values, then exit 
the program, your delay will be set to inf. 


References: 


Johnston, R.L., "NUMERICAL METHODS, a Software Approach", John 
Wiley and Sons, 1982 , 


Anton, Howard, “Elementary Linear Algebra", John Wiley and Sons, 
1981 


Atkinson, Kendal E., "An Introduction to Numerical Analysis", 
John Wiley and Sons, 1978 Z 
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USER INSTRUCTIONS 


Comments Input Display 
[------------+-+------------------ [-------------- [--------------------- I 
1) Run the program. Real or Complex? 

2) For real values: R 
For complex values: Cc order? 
3) Key in the number of rows: <value> A(1,1)=<value> 


4) 


5) 


6) 


The matrix must be square. 


You are now in the matrix 
editor. Please refer to that 
section, then continue with 
step 4. 


If you want to perform some 
operations and then make a 
few changes to the matrix, 
keep a copy of the original: 
If no copy is desired: 


Choose desired action. 

To create a new matrix: 
Cont. at step l. 

To edit the current matrix: 
Refer to the editor 
instructions, then cont. 
at step 4. 

To perform matrix operations: 
Continue at step 6. 

To quit the program: 


Calculation options. 

To solve system of equations: 
You are now in the matrix 
editor. Enter the values of 
the B vector for Ax=B. When 
you are done, you will see: 
Continue at solve section. 

To find the determinant: 
Continue at the determinant 
section. 

To find the inverse: 

Cont. at the inverse secton. 

To return to the main menu: 
Continue at step 5. 

To quit the program: 
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or 
R(1,1)=<value> 


keep a copy? 


Newmat,Edit,Calc,Quit 


Real or Complex? 
A(1,1)=<value> 

or 
R(1,1)=<value> 
Solv,Det,Inv,Main,Quit 


<blinking prompt> 


B(1)=<value> 


calculating... 


calculating... 


calculating.. 
Newmat,Edit,Calc,Quit 


<blinking display> 


USER INSTRUCTIONS 


Comments Input Display 
[------------------------------- I-------------- [--------------------- I 
kkKkKKK SOLVE kkkkkk 
$1) If a solution exists: <value x(1)> 
These are the values of [ENDLINE] <value x(2)> 
the result vector. [ENDLINE] 


Continue at step S2. 

If a solution does not 
exist, a message will be 
displayed: 


Return to step 6. 


$2) Once the last value has been 
displayed, you will see: 


$3) If you want to solve for 
a new B vector: Y 
Enter the values according 
to the edit instructions. 
Continue at step Sl. 
If you are done: N 
Return to step 6. 


kkkkKE DETERMINANT ***&*** 


Dl) The determinant will be 
dislayed: 


Press any key to continue. <any key> 
Go to step 6. 
kkkkketk INVERSE kkekekkk 
Il) If the inverse exists: 
Cont. at step I2. 
If no inverse exists: 
Return to step 6. 
I2) To list result by columns: Cc 
[ENDLINE] 
To list result by rows: R 
{ENDLINE] 


I3) Once all values have been 
listed, you will see: 
Go to step 6. 
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SINGULAR SYSTEM 
Newmat,Edit,Calc,Quit 


solve for new B? 


B(1) =<value> 


Solv,Det,Inv,Main,Quit 


det = <value> 
Solv,Det,Inv,Main,Quit 


list by Col or Row? 


SINGULAR SYSTEM 
Newmat,Edit,Calc,Quit 


<value A(1,1)> 
<value A(1,2)> 


<value A(1,1)> 
<value A(1,2)> 


Solv,Det,Inv,Main,Quit 


‘ 


USER INSTRUCTIONS 
kkekkkk MATRIX EDITOR kkkkkk 


The matrix editor allows the user to move through a matrix 
and change element values. For complex valued matrices, the 
real and complex parts of an element are edited one at a time. 


Movement through a matrix is accomplished with the arrow keys 
(left, right, up, and down), and element indices input. The 
movement wraps around when the boundary of a row or column is 
passed. 


Comments Input Display 
J -------- +--+ = - - = - - - - I--------~------ I --------------------- I 
El) To move from the current A(I,J) 
position: 

Move left: < A(I,J-1)= <value> 

Move right: > A(I,J+l)= <value> 

Move up: <up arrow key> A(I-1,J)= <value> 

Move down: <dwn arrow key> A(I+1,J)= <value> 


E2) To move to a desired 


element: {SPC] enter ROW,COLUMN 
<row>,<column> 
[ENDLINE] A(<row>,<col>)= <val> 
E3) To quit the editor: Q 
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EXAMPLES 


Example l. 


Find the determinant and inverse of the matrix below. It is 
assumed that the display is fix 4. 


6 3 -2 2 3 


1 4 -3 4 2 


Comments Input Display 
Te r  rn rrrnrrrrn I ------- = T -- enn ne I 
1) Run the program. , Real or Complex? 
2) The values are real. R order? 
3) Key in the number of rows: Sievoiwe] A(1,1)= 9.0000 
(avoiinwe] © 
6 
{ENDLINE] A(1,2)= 8.909 
(ENDLINE] ? 
3 
Enter the values of [ENDLINE] A(1,3)= 8.89800 
the matrix and use the [ENDLINE] ? 
editing features to adjust -2 
any incorrect values. ([ENDLINE] A(1,4)= 0.900 
{ENDLINE] ? 
2 
{ENDLINE] A(1,5)= 0.08909 
[ENDLINE] ? 
3 
{ENDLINE] A(5,4)= 0.090 
{ENDLINE] ? 
6 
[ENDLINE] A(5,5)= 8.9009 
[ENDLINE] ? 
2 
[ENDLINE] A(1,1)= 6.0009 
Q keep a copy? 
4) No copy will be needed. N Newmat,Edit,Calc,Quit 


5) Choose calculation 
options. Cc Solv,Det,Inv,Man,Quit 
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EXAMPLES 


calculating... 
det = -209.0009 


6) Calculate the determinant. D 


Press any key to continue. [ENDLINE]} Solv,Det,Inv,Main,Quit 
7) To find the inverse: I calculating.. 
list by Col. or Row? 
8) Choose to list by column. Cc C(1,l)= 9.2009 
{ENDLINE] C(2,1)=-1.7580 
[ENDLINE] C(3,1l)= 8.7988 
[ENDLINE] C(4,1)= 1.7258 
LENDLINE] C(5,1)= 1.9980 
[ENDLINE] C(1,2)=-0.1290 
[ENDLINE] C(2,2)=-2.9809 
{ENDLINE] C(3,2)= 1.2880 
[ENDLINE} C(4,2)= 2.5480 
[ENDLINE] C(5,2)= 1.4980 
{ENDLINE] C(1,3)=-9.9488 
{ENDLINE] C(2,3)= 0.5088 
{ENDLINE] C(3,3) =-G.2488 
[ENDLINE] C(4,3)=-0.5708 
[ENDLINE] C(5,3) =-8.2008 
[ENDLINE] C(1,4)= 6.9690 
[ENDLINE] C(2,4)= 1.7508 
4 {ENDLINE] C(3,4)=-0.5000 
[ENDLINE] C(4,4) =-1.6258 
[ENDLINE] C(5,4)=-1.0009 
{ENDLINE] C(1,5) =-3.3333E-12 
[ENDLINE] C(2,5)=-1.5088 
{ENDLINE] C(3,5) =-1.0908 
{ENDLINE] C(4,5)=-1.7589 
[ENDLINE] C(5,5)=-1.980G 
[ENDLINE] Solve,Det,Inv,Main,Quit 
9) Now find the determinant. D calculating.. 
det = -200.0000 
[ENDLINE] Solve,Det,Inv,Main,Quit 
18) Return to the main menu 
to perform the next problem. M Newmat, Edit,Calc,Quit 


Problem 2. 
Solve the following system of complex valued equations. 


2 + 3i 7 - li zl 2 + 21i 


4 - 1.3i 4+ Gi z2 1 + 3i 
Because the matrix created in problem 1 is real valued, a new 
matrix must be created. Therefore the example problem must start 


with the new matrix option. 
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1) 


2) 
3) 


4) 


5) 
6) 
7) 
8) 


9) 


EXAMPLES 


Comments 
wee ee I 


Start by creating a new 
matrix. 


Choose complex. 
The order is 2. 


You are now in the 
editor. 

The 'R' indicates that 
the real portion of 
element (1,1) is being 
prompted for. 

'I' indicates that the 
imaginary portion is 
being prompted for. 


Continue entering the 
values in the same manner. 
Refer to the editor 
instructions for editing 
features. 


When done, press: 

Keep a copy this time. 

Move to calculation options. 
Choose solve option. 

You are now in the editor.The 
prompt is for the real portion 


of element (1,1) from the 
b vector ( Ax = b). 


Continue until all 
values are correctly entered. 


18) To quit editing: 


You may use the arrow 
to view the entire display. 


2 
{ENDLINE] 


{ENDLINE] 
2 
[ENDLINE] 
[ENDLINE] 
3 
[ENDLINE] 
[ENDLINE] 
7 
[ENDLINE] 


e 
e 
e 
e 
° 


Q 
x 


{ENDLINE] 
2 
[ENDLINE] 
[ENDLINE] 
21 
(ENDLINE] 


Q 


[ENDLINE] 
[ENDLINE] 
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Display 


Real or Complex? 


order? 


R(1A)= 0.8908 


~~ 


I(1,1)= -6.9909 
? 
R(1,2)= 3.8008 
z 


I(1,2)= 2.0000 


keep a copy? 
Newmat,Edit,Calc,Quit 
Solv,Det,Inv,Main,Quit 


BR(1l )=4.0000 


? 


BI(1)= 6.9080 
? 


BR(2)= 8.9090 


calculating.. 

X(1)= 4.3565, 1.4203 i 
X(2)=-4.5681, @.7456 i 
solve for new B? 


11) Let's solve for a new b 
vector. 


12) Let's exit. 


EXAMPLES 
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Input Display 
[-------------- [--------------------- I 
Y BR(1)= 1.09909 
[ENDLINE] ? 
9 
({ENDLINE] BI(1)= 21.9000 
[ENDLINE] ? 
-.22 
{ENDLINE] BR(2)= 1.4794 
[ENDLINE] ? 
-3.5 
{ENDLINE] BI(2)= 3.0000 
1 
[ENDLINE] BR(1L)= 9.900 
Q calculating.. 
X(1)= 0.4800,-2.0302 i 
[ENDLINE] X(2)=-0.6952, 2.4362 i 
[ENDLINE] solve for new B? 
N Solv,Det,Inv,Main,Quit 
Q <blinking display> 


8019 
G928 
8930 
884d 
6058 
G86G 
8878 


8880 
0899 


8199 
9116 
81290 


9136 
8140 
81598 
9168 
0176 
8189 


G19 
8200 
8210 
0229 
8236 
8249 
8256 
0269 
8278 
G28 


G298 
0380 


8316 
8320 
0339 
0340 
8350 
8368 
637 
G38 


8390 
8408 
6418 
8429 
G43 
8446 
8459 
8469 


PROGRAM LISTING 


! rev 1.9 

DEF FNKS (D@$,K@S) ! Pind which key was hit 

DISP D@S 

"KEY': KS=KEYS 

IF NOT POS(K@S$,KS) THEN “KEY" 

FNKS=KS 

END DEF 

'n': T=POS("RC",FNKS ("Real or Complex?","RC") ) 

INPUT “order?";N ! Get the number of rows and columns. 
N=T*N @ M=N ! Set dimensions based on type. 

DIM A(M,N) ,S(M) 

R2=8 @ C2=G 

‘R's IF NOT C2 THEN GOTO 'El' 

FOR I=l1 TO M 

FOR J=l TO N 

A(I,J)=Al(I,J) ! Read in original matrix if copy was made 
NEXT J 

NEXT I 

‘El': CALL EDIT(A(,) ,M,N,T) 

R2=@ ! Indicate that the determinant has not been calculated. 
IF FNKS (“keep a copy?","YN")="y" THEN C2=1 ELSE C2=@ 
IF NOT C2 THEN GOTO 'M' 

DIM Al (M,N) 

FOR I=l TO M 

FOR J=l1 TON 

Al(I,J)=A(I,J) ! Save a copy if desired. 

NEXT J 

NEXT I 

'': GOTO FNKS("Newmat,Edit,Calc,Quit","NECQ") ! main menu 
'C'; GOTO FNKS("Solv,Det,Inv,Main,Quit","SDIMQ") 

"S's DIM B(M,T),X(M,1) ! solve invocation. 

‘S1l': CALL EDIT(B(,) ,M,T,@) 

DISP “calculating.." 

IF NOT R2 THEN CALL DECOMP(A(,),M,N,S()) @ R2=1 

IF S(M)=@ THEN DISP "SINGULAR SYSTEM" @ GOTO '‘'M' 

CALL SOLVE(A(,) -X(,) -B(,)/S() -N) 

CALL LIST(X(,) ,M,1,T,"X") 

IF FNKS("solve for new B?","YN")="Y" THEN "S1" ELSE "Cc" 
'r's DIM C(M,N) ,B(M,1),X(M,1) ! Inverse calculation 
DISP "calculating..." 

IF NOT R2 THEN CALL DECOMP(A(,),M,N,S()) @ R2=1 

IF S(M)=@ THEN DISP "SINGULAR SYSTEM" @ GOTO ‘'M' 

FOR J=1 TO N 

FOR I=l TO M 

B(I,1)=@ 

NEXT I 
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PROGRAM LISTING 


8478 B(J,1)=1 

9488 CALL SOLVE(A(,) ,X(,),B(,),S(),N) 
8499 FOR I=1 TO M 

8508 C(I,J)=X(I,1) 

@518 NEXT I 

0520 NEXT J 

G53@ CALL LIST(C(,),M,N,T,"C") 

8549 GOTO 'C' 


0550 ‘D's: IF T=2 THEN DISP "NOT DONE FOR COMPLEX" @ GOTO 'C' ! det. 
@568 DISP “calculating.." 

9570 IF NOT R2 THEN CALL DECOMP(A(,),M,N,S()) @ R2=1 

@58@ D=A(1,1) 

6590 FOR I=2 TOM 

0600 D=D*A(I,TI) 

0618 NEXT I 

@628 D=S(M)*D 

66380 DISP "“det= "sD 

@640 AS=KEYS @ IF AS="" THEN 640 ELSE "C" ! display till key hit 


650 'Q': PUT "#38" @ END ! restore blinking prompt and end. 


MATRIX EDITOR. Allows a matrix of dimension MxN to be edited. 

The type (T) can be 1 or 2. 1 indicates real, 2 indicates complex. 
If T=8, then the routine assumes a vector has been passed. The 
value of T is then changed to the correct type indicator value. 


0660 SUB EDIT(A(,) ,M,N,T) 

08670 DEF FNKS$(D@$,K@$) ! Get which key was hit 
0688 DISP D@S 

0698 "KEY': KS=KEYS 

@700 IF NOT POS(K@S,KS) THEN "KEY" 

0710 FNKS=KS 

0720 END DEF 


0738 DEF FNFS(Y) ! temporarily change display setting for index vals. 
8740 D9S=PEEKS ("2F6DC",2) @ STD 

0750 FNFS=STRS(Y) 

@76@ POKE "“2F6DC",D9S 

0776 END DEF 


@78@ DEF FNDS ! Create array elemnt prompt 

8790 IF NOT (T=1 OR MOD(J,2)) THEN DS=SS ELSE DS=RS 

0808 D1S=DS&" ("S&ENFS ( (I+T-1) /T) 

0810 IF DS[1,1]<>"B" THEN D1S=D1S&","&FNFS (INT ((J+T-1)/T) ) 
0828 IF T=2 AND NOT MOD(J,2) THEN V=-A(I,J) ELSE V=A(I,J) 
0830 FNDS=D1S&")= "&STRS(V) 

0849 END DEF 


8858 I=1 @ J=1 @ R=2 @ TS="BRBI" 

0860 IF T=2 THEN RS="R" @ SS="I" 

0870 IF T=1 THEN RS="A" 

6880 IF T=0 THEN RS=TS[1,N] @ SS=TS[3,N+2] @ T=N 
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8898 


8989 
6919 
8928 


8936 
8946 


PROGRAM LISTING 


Figure out which key has been hit and branch to legal choice. 
'LOP': AS=FNKS(FNDS,"Q #38#47#48#50#51") @ 

IF POS ("9134578#",AS) THEN "LOP" 
IF LEN(AS$)=1 THEN GOSUB UPRCS (CHRS (NUM(AS$)+33)) @ GOTO "LOP" 
AS {L Ly="K" 
GOSUB AS @ GOTO "LOP" 


'K38's INPUT A(I,J) @ IF T=1 THEN GOTO 'K48' ! enter a value 
IF MOD(J,2) THEN A(I+1,J+1)=A(I,J) ELSE A(I+l,J-1)=A(I,J) @ A(I, 


J)=-A(I,J) 


8956 
8966 
G97 


8988 
8996 
1606 


1916 
1926 
1830 


1948 
1958 
1868 


1076 
1986 
1696 
1168 
11196 
1128 


11390 


'k48': J=J+l @ IF J<=N THEN RETURN ! move right 
J=1 @ I=I+T @ IF I>M-T+l THEN I=1 
RETURN 


'K47': J=J-1 @ IF J THEN RETURN ! move left 
J=N @ I=I-T @ IF I<l THEN I=M-T+1 
RETURN 


"K5@': I=I-T @ IF I>@ THEN RETURN ! move up 
I=M-T+l1 @ J=J-1 @ IF NOT J THEN J=N 
RETURN 


"K51': [=I+T @ IF I<=M-T+1 THEN RETURN ! move down 
I=l @ J=J+1l @ IF JON THEN J=1 
RETURN 


‘A's ON ERROR GOTO 'A' ! user specified move 
INPUT "enter ROW,COLUMN";I,J 

IF I<l OR J<l THEN I=M+1 

IF T-l1 THEN I=2*I-1 @ J=2*J-1 

IF I>M OR J>N THEN DISP "OUT OF BOUNDS" @ GOTO 'A' 
OFF ERROR @ RETURN 


"R': POP @ SUB END 


LIST MATRIX SUB PROGRAM. Allows a matrix to be listed by row 
or column. The array name is passed through parameter BS. The 
matrix is MxN, and the type (T) is 1 for real values, and 2 


for complex values. If a vector is passed, the routine will 
list by column. A vector is implied by BS = 'X'. 
1148 SUB LIST(A(,) ,M,N,T,BS) 
115@ DEF FNFS(Y) ! Create index integer prompts 
1168 D9S=PEEKS("2F6DC",2) @ STD 
1170 FNFS=STRS (Y) 
1188 POKE "2F6DC",D9S 
1198 END DEF 
12860 DEF FND$S ! Create element prompt 
1210 IF BS="X" THEN FNDS=FNFS((I+T-1)/T) ELSE 
FNDS=FNFS ( (I1+T-1)/T) &","&FNFS ((J+T-1) /T) 
1228 END DEF 
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PROGRAM LISTING 


1230 D1S=PEEKS$("2F946",4) @ DELAY INF,INF ! Temporarily change delay 
1240 IF BS="X" THEN GOTO 'C' 


1258 'P': DISP "list by Col. or Row?" ! get display choice 
1260 AS=KEYS$ 

1278 IF NOT POS("CR",UPRCS(A$[1,1])) THEN 'P' 

1286 GOTO UPRCS(AS[1,1]) 


1299 'C': FOR J=l TO N STEP T ! display by column 
130@ FOR I=l1 TO M STEP T 

1318 DISP BS&" ("&FNDS&")="ZA(I,I) 5 

1320 IF T=2 THEN DISP ",";A(I+1,J) ;"i"; 

1338 DISP 

1348 NEXT I- 

1358 NEXT J 

1368 I=l @ J=1 @ GOTO 'E' 


1370 'R': FOR I=l TO M STEP T ! display by row 
1388 FOR J=l1 TO N STEP T 

1398 DISP BS&"("&FNDS&")="ZA(I,J); 

14060 IF T=2 THEN DISP ","sA(It+1,J);"i"; 

141@ DISP 

1428 NEXT J- 

1438 NEXT I 

1440 I=1 @ J=1 


145@ 'E': POKE "2F946",D1$ @ SUB END ! restore delay and exit 


DECOMPOSITION OF MATRIX. Performs an LU decomposition of an MXN 
matrix using partial pivoting. The pivoting strategy is recorded 
in vector S. 


1466 SUB DECOMP(A(,) ,M,N,S()) 
1478 S(M)=1 

1488 FOR R@=1 TO M-1 

1498 PO=RD 

150@ P1l=A(R@,RG) 


1510 FOR I=R@+1 TO M ! choose largest absolute value for pivot 
1520 IF ABS(A(I,R@))>ABS(P1) THEN P@=I @ P1=A(I,RG) 

1530 NEXT I 

1548 IF A(PG,R@)=@0 THEN S(M)=8 @ GOTO 'END' ! quit if singular 
1558 S(RB)=Pd 

1560 IF P@=R8 THEN GOTO 'C' 


1578 FOR I=R@ TO N ! row exchange 
1589 T=A(RQG,I) 

159@ A(R@,1I)=A(PG,I) 

1608 A(PG@,1I)=T 

16198 NEXT I 

1628 S(M)=-S (M) 


163@ 'C': FOR R1=RG+1 TO M ! row rl <--rl-mult*rg 
4 


page 41 


PROGRAM LISTING 


1640 M1L=A(R1,RG0)/A(RGO,RG) ! form multiplier ( 
1658 A(R1,RO)=Ml ! and save it. 

1660 FOR E=R9+1 TO N 

1678 A(R1,E) =A(R1,E) -M1*A(RG@,E) 

1688 NEXT E 

1698 NEXT R1 

1708 NEXT RG 

1710 IF A(M,N)=8 THEN S(M)=@ 

1728 ‘END': SUB END 


SOLVE ROUTINE. Takes an LU form matrix, pivot strategy vector, B 
vector, and calculates the X vector for the matrix equation 
Ax=b. 


173@ SUB SOLVE(A(,),X(,),B(,),S() ,N) 
1740 M=N 


1750 FOR I=1 TO M-1 ! Permute B and perform reduction 

176@ T=B(S(I),1) 

1770 B(S(I),1)=B(I,1) 

1788 B(I,1)=T 

1799 FOR J=I TO M-1 

1800 B(J+1,1)=B(J+1,1)-B(1I,1) *A(J+1,1) 

1810 NEXT J : 
1828 NEXT I ( 


1838 FOR I=N TO 1 STEP -1 ! back substitution 
1848 X(I,1)=B(I,1) 

1858 FOR J=I+l TO M 

1868 X(I,1)=X(1I,1)-A(1I,J) *X(J,1) 

1870 NEXT J 

188@ X(I,1)=X(1I,1)/A(I,TI) 

1898 NEXT I 

1998 SUB END 


page 42 


Fourier Transforms 


—————————————————— —— — — — —— —  ——eeeeeeeeeeeeeeeeEOoeoeoeoeoeoeoeoeoeoeq®q®o _ee_eeeeeeeess eee SSS SS S——ww_as 


This program calculates a fast Fourier transform from a set of time 
domain points to a set of frequency domain points. The inverse fast 
Fourier transform, calculating the set of time domain points from a 
set of frequency domain points, may also be calculated. 


The method used is a modification of the basic FFT algorithm. The 
modified algorithm takes advantage of the fact that series data is 
real, and uses the space normally reserved for the imaginary part of 
the complex sequence to calculate a double-length real transform. 
This is represented for two "N" length transforms as: 


Z(n) = X(n) + iY(n) @<n < N data points 


The transform is: 


Z(m) = X(m) + iY¥(m) 
Z(m) + Z(N-m)* 
where X(m) = ---------------- 
2 
Z(m) - Z(N-m) * 
Y(m) = ---------------- 


zZ* is the complex conjugate of Z. 
The time series F(n) is given by: 


F(n) = X(2n) + Y(2n+1) 


The transformation of this is: 


N-1 N-1l 
F(m) = > X(2n)w*mn + Ps Y (2n+1)w*mn 
n= n=@ 
N-1 N-l 
= SO X(p)w*2mp + Y ¥(p) (w°2mp) (w*m) 
p=9 p=9 
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and: 


F (m) X(m) + Y¥(m)w*m (1) 


F (N-m) X*(m) - [¥(m)w*m]* (2) 


Similarly, the inverse transform may be obtained from (1) and (2): 


F(m) + F(N-m)* F(m) -— F(N-m) * 
Z(mM) = ---------------- + iw (-m) ---------------- 
2 2 


F(m) + F(N-m)* * F(m) - F(N-m)* * 


This is simply an interchange of Z(m) and F(m) in (1) and (2), and 
substitution of -w*(-m) for wm. 


The advantages gained from this adaptation of the general FFT 
algorithm for time series data are: 


(a) A transform of twice the length can be handled with no increase 
in storage for input data. 


(b) Since the calculation of the transform is structured as an 
interactive process, intermediate and final results are stored in 
the same locations used for input. 

NOTE: 


1. Since F(@) and F(N) are real only, F(N) can be stored in the 
imaginary location of F(@), i.e., £f(1). 


2. wm = c*(-2im*pi/2N). This is half the minimum value of rotation 
normally used in an N-point transfer. 


3. * denotes the complex conjugate. 


REFERENCES 
Brigham, E. O., THE FAST FOURIER TRANSFORM, Prentice-Hall, Inc. 1974. 


FAST FOURIER TRANSFORM, HP-85 Math Pac, Hewlett-Packard, 1979. 
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USER INSTRUCTIONS 


Comments Input Display 
[---------------------- [-------------- [--------------------- I 


1) Run the program. TIME/FREQ. DATA? (T/F) 


2) If time data is to be input, 
press T and continue with 
step 3. If frequency data is 
to be input, go to step ll. {[T] #* OF DATA POINTS? 


***TIME DOMAIN DATA*** 


3) Enter number of time domain 
data points to be used. Must 
be an integer power of 2 and 
>2. If it is not, the message 
"INPUT OUT OF RANGE" will be 
displayed and you will be 
prompted to enter the number 
again. If available memory is 
insufficient for the specified 
# of points, the message "NOT 
ENOUGH MEMORY" will be dis- 
played and the program will 
again prompt for the number 
of points. <N> [ENDLINE] DATA POINT (nnn) ? 


4) The program will now prompt 
for data points 1 through N, 
where N is the # specified 
in 3). <value>[ENDLN] DATA POINT (nnn) ? 
<value>[ENDLN] . 


. DATA POINT( N )? 
<value> [ENDLN] 


5) If any mistakes were made in CHANGES? (Y/N) 
input, you now have the 
opportunity to correct the 
data. Press [Y] to make any 
changes, or [N] to go to 


step 6). {Y] DATA POINT TO CHANGE? 
Enter point # to change. <nnn> 
[ENDLINE] <old value> 
Original entry will be 
displayed for reference. <new value> 
[ENDLINE] CHANGES? (Y/N) 


If no (further) changes are 
necessary, press [N]. (N] TRANSFORMING... 
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USER INSTRUCTIONS 


Comments Input Display 
[-----~------------------------- [-------------- [--------------------- I 


6) Calculation of FFT has begun. DC TERM = 


7) Calculation ends and DC term 
is displayed. Press any key 
to continue. <any key> MAX FREQ. = 


8) Display maximum frequency. 
Press any key to resume 
output. <any key> FREQ. DOMAIN OUTPUT: 


9) Display complex data pairs of 

calculated frequency domain. 

First the real, and then the 

imaginary part of the data 

pair will be displayed. Press 

any key to display each value. LR=<value> 
<any key> 1lI=<value> 
<any key> 2R=<value> 
<any key> . 


. NR=<value> 
<any key> NI=<value> 
<any key> > 

18) For a new problem, go to 

step l. 


***FREQUENCY DOMAIN DATA*** 


11) If frequency data is to be 
input, press [F]. [F] # OF COEFF. PAIRS? 


12) Enter the number of coeffic- 
ient pairs. This number must 
be one less than an integer 
power of 2 (1,3,7,...). If 
it is not, the message 
“INPUT OUT OF RANGE" will be 
displayed and you will be 
prompted to enter the number 
again. If available memory is 
insufficient for the number 
of pairs specified, the 
message "NOT ENOUGH MEMORY" 
will be displayed and the 
program will again ask for 
the number of coefficient 
pairs. <N> [ENDLINE] DC TERM= 


13) Enter the DC term. <DC term> MAX FREQ. TERM= 
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USER INSTRUCTIONS 


Comments Display 


14) Enter the maximum frequency 


term. <max freq> FREQ. DOMAIN DATA - 
15) Enter the real and imaginary 
coefficients for each data REAL( 1 )? 
pair as prompted. <real coeff.> IMAG( 1 )? 
<imag coeff.> REAL( 2 )? 
<real coeff.> ; 
° REAL( N )? 
<real coeff.> IMAG( N )? 
<imag coeff.> CHANGES? (Y/N) 


16) If any mistakes were made 
in entering the coefficient 
pairs, they may be corrected 
by pressing [Y]. If no changes 
are necessary, press [N] and 


go to step 17). DATA PAIR TO CHANGE? 


[Y] 
Enter # of coefficient pair 


to correct. REAL( n )<old entry> 


<n> [ENDLINE] 


<new value> 
[ENDLINE] 


Original entry is displayed 


for reference. IMAG( n )<old entry> 


<new value> 


{ENDLINE]) CHANGES? (Y/N) 
If no (further) changes are 
necessary, press [N]. (N] TRANSFORMING... 
17) Calculation of inverse FFT 
has begun. TIME DOMAIN OUTPUT: 
18) Calculation ends and data 
points are displayed. Press 
any key to move to next 
data point. PT( 1 )=<value> 
<any key> PT( 2 )=<value> 
<any key> ‘ 


a PT( N )=<value> 
<any key> > 


19) For a new problem, 
step 1). 


go to 
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EXAMPLES 


A) For the following set of time domain data points, calculate the 
Fourier transform to frequency data. 


P(1)=1 P(9)=-1 
P(2)=1.3066 P (10) =-1. 3066 
P(3)=1.4142 P(11)=-1.4142 
P(4)=1.3066 P (12) =-1. 3066 
P(5)=1 P(13)=-1 
P(6)=.5412 P(14)=-.5412 
-P(7)=0 P(15)=9 
P(8)=-.5412 P(16)=.5412 


USER INSTRUCTIONS 


Comments Input Display 
J ee a nr nn rrr rrr rrr J -------------- J ----- - ee eee I 
1) Run the program. TIME/FREQ. DATA? 
2) Choose time data input. {T] # OF DATA POINTS? 
3) Enter # of points. 16 [ENDLINE] DATA POINT( 1 )? 
4) Enter point values. 1 [ENDLN] DATA POINT( 2 )? 
1.3066 [ENDLN] DATA POINT( 3 )? 
1.4142 [ENDLN] DATA POINT( 4 )? 
1.3866 [ENDLN] DATA POINT( 5 )? 
1 [ENDLN] DATA POINT( 6 )? 
25412 [ENDLN] DATA POINT( 7 )? 
g [ENDLN] DATA POINT( 8 )? 
-.5412 [ENDLN] DATA POINT( 9 )? 
-1 [ENDLN] DATA POINT( 19 )? 
-1.3066 [ENDLN] DATA POINT( ll )? 
-1.4142 [ENDLN] DATA POINT( 12 )? 
-1.3066 [ENDLN] DATA POINT( 13 )? 
-1 [ENDLN] DATA POINT( 14 )? 
-.5412 [ENDLN] DATA POINT( 15 )? 
G [ENDLN] DATA POINT( 16 )? 
-5412 [ENDLN] 
5) FFT calculation begins. TRANSFORMING... 
6) Frequency domain output. FREQ DOMAIN OUTPUT: 
DC TERM=9 
<any key> MAX FREQ. =@ 
<any key> FREQ Goria ort Pett 
<any key> 1R=-1.000010E+000 
<any key> 2R= GO. GGGIBDE+OOG 


<any key> 2I= 8. GBOOOE+OGG 
<any key> 3R=-1.339454E-006 
<any key> 31=-1.339454E-806 
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Comments Input Display 
J on rrr rrr [ween ee ---- -- [------- ee I 
<any key> 4R= G8. GGBGGTE+GO9 
<any key> 4Il= G8. 9DPBOBE+BO0 
<any key> 5R= 6.134477E-806 
<any key> 51=-6.134477E-906 
<any key> 6R= GB. GGGGDGE+O9G 
<any key> 6I1= G6. GGDGBOE+BGD 
<any key> 7R=-1.592234E-0885 
<any key> 71=-1.502234E-965 
Done. <any key> > 


B) 


Comments Input Display 
[------------------------------- I-------------- [--------------------- I 
1) Run the program. TIME/FREQ. DATA? (T/F) 
2) Enter frequency data. {F] # OF COEFF. PAIRS? 

3) Seven frequency data pairs. 7 ([ENDLINE] DC TERM= 
4) DC term is @. @  [ENDLINE] MAX FREQ. TERM= 
5) Maximum frequency term is @. @  [ENDLINE] FREQ. DOMAIN DATA - 
6) Begin entry of frequency REAL( 1 )? 
domain data pairs. 1 [ENDLINE] IMAG( 1 )? 
-1 [ENDLINE] REAL( 2 )? 
G [ENDLINE] IMAG( 2 )? 
Gg [ENDLINE] REAL( 3 )? 
-1.3395E-6 
[ENDLINE] IMAG( 3 )? 
-1.3395E-6 
(ENDLINE] REAL( 4 )? 
g (ENDLINE] IMAG( 4 )? 


USER INSTRUCTIONS 


For the following set of frequency domain data pairs, perform the 
inverse Fourier transform to calculate the set of time domain data 
points. The DC term and the maximum frequency term are @. 


REAL (1) =1 
REAL (2) =@ 


REAL (3) =-1.3395E-6 


REAL (4) =@ 


REAL (5) =6.1345E-6 


REAL (6) =9 


REAL (7) =-1.5022E-5 


IMAG (1) =-1 


IMAG (2) =@ 
IMAG (3) =-1.3395E-6 
IMAG (4) =@ 
IMAG (5) =-6.1345E-6 
IMAG (6) =@ 
IMAG (7) =-1.5@822E-5 


USER INSTRUCTIONS 
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7) 


8) 


Inverse 
begins. 


Display 
points. 


Done. 


FFT calculation 


time domain data 


INSTRUCTIONS 


Input Display 
[-------------- [--------------------- I 
g {ENDLINE] REAL( 5 )? 
6.1345E-6 
{ENDLINE] IMAG( 5 )? 
-6.1345E-6 
[ENDLINE] REAL( 6 )? 
G (ENDLINE]} IMAG( 6 )? 
G {ENDLINE] REAL( 7 )? 
-1.5822E-5 
{ENDLINE] IMAG( 7 )? 
-1.59822E-5 
{ENDLINE] 
TRANSFORMING... 
TIME DOMAIN OUTPUT: 
PT( 1)= 9.999898E-901 
<any key> PT( 2)= 1.306587E+00G 
<any key> PT( 3)= 1.414186E+009 
<any key> PT( 4)= 1.306587E+008 
<any key> PT( 5)= 9.999898E-601 
<any key> PT( 6)= 5.411945E-@01 
<any key> PT( 7)= 1.000000E-G12 
<any key> PT( 8) =-5.411945E-001 
<any key> PT( 9)=-9.999898E-801 
<any key> PT( 10) =-1.306587E+009 
<any key> PT( 11)=-1.414186E+080 
<any key> PT( 12) =-1.306587E+000 
<any key> PT( 13)=-9.999898E-081 
<any key> PT( 14) =-5.411945E-061 
<any key> PT( 15)=-1.90000GE-G12 
<any key> PT( 16)= 5.411945E-@01 
<any key> > 
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G91 
8920 
8838 
GB4B 
8058 
G86G 
8978 
8880 
8898 
G10 
G118 
8120 


keke 
kkk 
kkkKK 


kkk 


9139 
9140 
91590 
81608 
6178 
0188 
8196 
820 
8216 
0220 
0238 
G249 
8258 
8268 
6278 


kkkek 


LISTING 


! FOURIER TRANSFORM 
! Revision 1.0 4/16/84 
! 


F9=FLAG (5,FLAG(-19) ) 

OPTION BASE 1 @ OPTION ANGLE RADIANS @ STD @ DELAY 9,0 
ON ERROR GOTO 'ERR' 

I1$="TF" @ 12$="YN" 

O1$="3D,'R=',MZ.6DE" @ 02$="3D,'I=",MZ.6DE" 
03S="'"PT(',3D,') =" ,MZ.6DE" 

DISP 'TIME/FREQ. DATA? (T/F)' 

'WL': KIS=UPRCS$(KEY$) @ IF NOT POS(I1$,K1$) THEN ‘W1' 
IF K1$='F' THEN SFLAG 1 @ F=-1 ELSE CFLAG 1 @ Fel 


Input data points for FFT, or coefficient pairs for **** 
inverse FFT. Number of data points must be a power eeRE 
of 2 and greater than 2. Number of coefficient weaK 
pairs must be one less than a power of 2. eRe 


"IN's 

IF NOT FLAG(1) THEN INPUT '# OF DATA POINTS?';N @ GOTO ‘INA' 
INPUT '# OF COEFF. PAIRS?';N 

N=N*¥2+2 

"INA'Ts: Pe=l 

IF N=2 THEN 'OR' 

FOR L=l1 TO 1@ 

P=P*2 

IF P=N THEN Pl=L @ N2=N/2 @ GOTO 'START' 

NEXT L 

‘OR': 

DISP 'INPUT OUT OF RANGE' @ WAIT 1 @ GOTO 'IN' 
'"START': 

DIM R(N2),I(N2) 

IF FLAG(1) THEN ‘IFFT! 


Input time domain data for FFT **** 


'FFT': 

J=G 

FOR L=l1 TO N2 

J=J+1l @ DISP 'DATA POINT(';J;')'; @ INPUT R(L) 
J=J+1 @ DISP 'DATA POINT(';J;')'; @ INPUT I (L) 
NEXT L 


Changes to FFT data **** 


"Crs 

DISP 'CHANGES? (Y/N)' 

'w2': K2S=UPRCS(KEYS) @ IF NOT POS(I2S,K2$) THEN 'W2' 
IF K2S="N" THEN DISP 'TRANSFORMING...' @ GOTO ‘FFTC' 
'C2's INPUT 'DATA POINT TO CHANGE?';L 

IF L<=@ OR L>2*N2 THEN 'C2' 

I2=L/2 

IF I2=INT(12) THEN J=1I2 ELSE J=INT(I2)+1 

DISP "DATA POINT(';L;')'; 
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8439 
G44 


kkekk 
kkkk 


8458 
8460 
8478 
8480 
8498 
8590 
85190 
8520 


kkRK 


8536 
0540 
0558 
8560 
8570 
8580 
8598 


LISTING 


IF J=1I12 THEN INPUT '',STRS(I(J));1(J) ELSE INPUT ++ sTRS(R(J)) RIT 
GoTo 'Cl' 


Input DC term, maximum frequency and eaKK 
frequency domain data for inverse FFT. **** 


‘IFFT': 

INPUT 'DC TERM=';R(1) 

INPUT 'MAX FREQ. TERM=';I (1) 

DISP ‘FREQ. DOMAIN DATA -' @ WAIT 1 


FOR L=2 TO N2 

DISP 'REAL(';L-1;')'; @ INPUT R(L) 
DISP 'IMAG(';L-1;')'; @ INPUT I(L) 
NEXT L 


Changes to inverse FFT data **** 


'C3': 

DISP 'CHANGES? (Y/N)' 

'W3': K3S=UPRCS(KEYS) @ IF NOT POS(I2$,K3$) THEN 'W3' 
IF K3$='N' THEN DISP 'TRANSFORMING...' @ GOTO 'IFFTC' 
'c4's: INPUT 'COEFF. PAIR TO CHANGE?';L 

IF L<=@8 OR L>N2-1 THEN *C4' 

DISP 'REAL(';L;")'; @ INPUT '',STRS(R(L+1)) ;R(L+1) 


0696 DISP 'IMAG(';L;')'; @ INPUT '',STRS(I(L+1));1I(L+1) 

9619 GOTO 'C3' ( 
keke «Start FFT calculation **** 

6620 ‘FFTC': 

0630 K=6 

064% FOR J=1 TO N2-1 

0650 L=2 

@66@ IF K<N2/L THEN 689 

@670 K=K-N2/L @ L=L+L @ GOTO 669 

@680 K=K+N2/L 

@690 IF K<=J THEN 719 

0700 A=R(J+1) @ R(J+1)=R(K+1) @ R(K+1)=A @ A=I(J+1) @ 1(J+1)=1(K+1) 
@ I(K+1)=A 

6718 NEXT J 

0720 G=.5 @ P2=1 

8730 FOR L=l TO Pl-1l 

0740 G=G+G @ C=1 @ E=6 @ Q=SQR((1-P2)/2)*F 

0750 P2=(1-2* (L=1) ) *SQR((1+P2) /2) 

0760 FOR M=l1 TO G 


8776 
8786 
8796 
G80G 
8818 
8820 
8830 
8848 


FOR J=M TO N2 STEP G+G 
K=J+G @ A=C*R(K)+E*I(K) @ B=E*R(K)-C*I (K) 


R(K)=R(J)-A @ I(K)=I1(J)+B @ R(J)=R(J)+A @ I(J)=1(J)-B 

NEXT J 

A=E*P2+C*Q @ C=C*P2-E*Q @ E=A 

NEXT M 
NEXT L 
IF FLAG(1) THEN ‘IFFTO' 
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kak 


G85 
8869 
8878 
G88 
8899 
6980 
89198 
8928 
8936 
8949 
89598 
8969 
8970 
8980 
8998 
19090 


kkk 


1919 
102 
18398 
1049 
1958 
1869 
1876 
198@ 
1996 


keKK 


1196 
1116 
11296 
1136 
1140 
1159 
1166 
1179 
1186 
1196 
1299 
12190 
1226 


LISTING 
Start inverse FFT calculation **** 


"TFFTC!: 

A=PI/N2 @ P2=COS(A) @ Q=F*SIN(A) 

A=R(1) @ R(1L)=AtI(1) @ 1(1)=A-I(1) 

IF NOT FLAG(1) THEN R(1)=R(1)/2 @ I(1)=1(1)/2 

C=F @ E=@ 

FOR J=2 TO N2/2 

A=E*P2+C*Q @ C=C*P2-E*Q @ E=A @ K=N2-J3+2 @ A=R(J)+R(K) 
B=(I(J) +1 (K))*C-(R(J)-R(K))*E @ U=I(J)-I(K) 
V=(1I(J)+1I(K)) *E+(R(J)-R(K) ) *C 

R(J)=(A+B)/2 @ I(J)=(U-V)/2 @ R(K)=(A-B)/2 @ I(K) =-(U+V)/2 
NEXT J 

I (N2/2+1) =-I (N2/2+1) 

IF FLAG(1) THEN: 'FFTC' 

FOR J=l1 TO N2 

R(J)=R(J)/N2 @ I(J)=1 (J) /N2 

NEXT J 


FFT output **** 


‘FFTO': 

DISP 'DC TERM =';R(1) @ GOSUB 'WAIT' 
DISP 'MAX FREQ. =';1I(1) @ GOSUB ‘WAIT’ 
DISP 'FREQ DOMAIN OUTPUT:' @ WAIT 1 
FOR L=2 TO N2 


DISP USING 01$;L-1,R(L) @ GOSUB ‘WAIT’ 
DISP USING 02$;L-1,I1(L) @ GOSUB ‘WAIT’ 
NEXT L 
GOTO 'DONE' 
Inverse FFT output **** 
"TFFTO': 
DISP 'TIME DOMAIN OUTPUT:' @. WAIT 1 
J=1 
FOR L=l TO N2 
J=J+1 
DISP USING 03$;J-1,R(L) @ GOSUB ‘WAIT' 


DISP USING 03$;J,I(L) @ GOSUB ‘WAIT’ 

J=J+1l 

NEXT L 

'DONE': F9=FLAG(-19,FLAG(5)) @ PUT '#43' @ END 

'WAIT': IF KEYS='' THEN 'WAIT' ELSE RETURN 

'ERR': IF ERRL=26@ THEN DISP 'NOT ENOUGH MEMORY' @ GOTO ‘IN! 
DISP ERRMS @ GOTO 'DONE' 
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Polynomial Root Finder 


This program finds all solutions, both real and COM N SR of P(x)=@, 
where P is a polynomial of the form: 








P(x) =a(n)x*n + a(n-)x*(n-1) + ... + a(1)x + a(@) = @ 


Inputs to the program are the degree of the polynomial, the real 
coeficients a(n)...a(@), tolerances for the evaluation of the function 
and for each root, and the maximum number of iterations per root. 


This program uses Laguerre's method to find the roots of the specified 
polynomial by computing a sequence of approximations Z(1l),_ 22) peooy 
to a root using the formula Z(k+1)=Z(k)+S(k). S(k) is called the 
Laguerre step, and is defined as: 


-nP (Z(k) ) 


P'(Z(k))+[(n-1) *2(P' (Z(k)) *2-n(n-1)P'' (z(k))1]7.5 


where P, P', and P'' are the value of the polynomial and its first and 
second derivatives evaluated at the current iterate k, and n is the 
degree of the polynomial. The sign in the denominator is chosen to 
give the Laguerre step of smaller size, which in most cases insures 
that the roots will be found in order of increasing magnitude. 


After an iterate is accepted as a root, synthetic division is used to 
deflate the polynomial by the factor (x-r) if the root is real, or 
(x*2-2Re(r)+!irt!*2) if the root is complex. This saves arithmetic 
operations, and prevents repetitive convergence to the same root. 


For polynomials with only real roots, Laguerre's method will always 
converge to a root for any choice of real initial estimate. However, 
for roots of high multiplicity, some loss of accuracy may be observed. 
If complex roots are present, this method will usually converge to a 
valid root. If it does not, provisions are made for supplying a new 
initial estimate and starting the process again. 

REFERENCES 

Dahlquist, G. and Bjorck, A. NUMERICAL METHODS. Prentice-Hall, 1974. 


HP-75 MATH PAC. Hewlett-Packard, 1983 
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USER INSTRUCTIONS 
COMMENTS INPUT ; DISPLAY 
[------------------------------- [-------------- [--------------------- I 
1) Run the program. POLYNOMIAL ROOT FINDER 
ORDER OF POLYNOMIAL? 
2) Input degree of polynomial. 
Must be a positive integer 
greater than 1. If it is not, 
you will be asked to enter 
it again. <n> A(n) =? 
3) Enter coefficients of each 
term, starting with the 
highest-ordered term. <a(n)> 
[ENDLINE] A(n-1)=? 
<a(n-1)> 
{ENDLINE] - 
: A(@) =? 
<a(@)> 
[ENDLINE] TOL. FOR ROOTS=1.E-18 
4) Input tolerance for roots. 
Default value of 1E-19 is 
displayed. If the magnitude 
of the Laguerre step is less 
that this value (and 5) is 
also satisfied) then the 
current iterate is accepted 
as a root. <new value> 
[ENDLINE] TOL. FOR FCN=1.E-8 
5) Input tolerance for evaluation 
of function. default value of 
1E-8 is displayed. If P(x) for 
the current iterate is less 
than this value, and step 4) 
has been satisfied, the current 
iterate is accepted as a root. <new value> 
{ENDLINE] MAX # OF ITERATIONS= 
6) Input maximum number of 
iterations for each root. If 
this number is exceeded before 
a valid root is found, the 
message 'NO CONVERGENCE' will 
be displayed and the user will 
be allowed to specify a new 
initial iterate and start the 
search again. <nnn> 
{ENDLINE] LOOKING FOR ROOTS... 


page 55 


USER INSTRUCTIONS 


COMMENTS INPUT DISPLAY 
[------------------------------- I-------------- [--------------------- I 
7) Calculation of roots has begun. # OF ROOTS FOUND = 1 

As each root is found, the # OF ROOTS FOUND = 2 
program will indicate how many ° 
roots have been found to this : 
point. ‘: 
# OF ROOTS FOUND = n 


ROOT# 1: R=<value> 
8) The real and imaginary parts 
of each root are displayed. 
Pressing any key will continue 


the displaying of the roots. <any key> ROOT# 1: I=<value> 
<any key> ROOT# 2: R=<value> 
<any key> ° 


. ROOT# N: R=<value> 
<any key> ROOT# N: I=<value> 
<any key> > 

9) Done. 
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EXAMPLES 


A) Find the roots of the polynomial given below. Use default values 
for the tolerances and limit iterations to 16. 


P(x) = 5x°6-45x°54+225x~4-425x"34+170x"2+370x-508 


USER 


COMMENTS 


INSTRUCTIONS 


DISPLAY 


[----+=.0-------- 5+ ------4--- [-----s5-5<---- [--------------------- I 


1) Run the program. 


2) Enter order of polynomial. 
3) Enter the coefficients, 


starting with the highest 
order term. 


4) Use default value. 
5) Use default value. 
6) Limit to 1@ iterations. 


7) Calculation of roots begins. 


8) Display real and imaginary 
parts of each root. 


9) Done. 


6 [ENDLINE] 


s [ENDLINE] 
-45 [ENDLINE] 
225 [ENDLINE] 
-@25 [ENDLINE] 
tr [ENDLINE] 
478 [ENDLINE] 
-50@ [ENDLINE] 


[ENDLINE] 
[ENDLINE] 


18 [ENDLINE] 


<any key> 
<any key> 
<any key> 
<any key> 
<any key> 
<any key> 


<any key>. 


<any key> 
<any key> 
<any key> 
<any key> 
<any key> 
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POLYNOMIAL ROOT FINDER 
ORDER OF POLYNOMIAL? 


# OF ROOTS FOUND 
# OF ROOTS FOUND 
# OF ROOTS FOUND 
# OF ROOTS FOUND 


ROOT # 


ROOT# 
ROOT# 
ROOT# 
ROOT# 
ROOT# 
ROOT# 
ROOT# 
ROOT# 
ROOT# 
ROOT# 
ROOT# 
> 


A(6) =? 


A(5)=? 
A(4) =? 
A(3)=? 
A(2)=? 
A(1)=? 
A(G) =? 
TOL. 


FOR ROOTS=1.E-19 


TOL. FOR FCN=1.E-8 
MAX # OF ITERATIONS= 


LOOKING FOR ROOTS... 


new 


R= 1. 000G0GE+800 


I= |. OBOGBGE+OGD 
= 1. 000000E+000 
=~-1.900DBGE+BGG 
=-]1.G9G0GSE+9G0 
= 2. 000000E+000 
= 2. 000000E+000 

I= 0. 0000G0E+R00 

R= 3.000000E+@00 
=-4.G0000GE+00O 

R= 3.900000E+900 

I= 4, 00000GE+009 


8910 
6828 
8938 
8B4G 
8059 
8860 
0878 
898 


kkk 
kkk 


8998 
81808 
81190 
8120 
8138 
8149 
8159 
8160 
81798 
8180 
81998 
82008 
6218 
8228 
82390 
8249 
8250 
8260 


kkk 


LISTING 


DEFAULT ON 

! POLYNOMIAL ROOT FINDER 

! Revision 1.00 4/8/84 

STD @ OPTION BASE @ @ INTEGER I,J,K,N 
A@S="'YN' @ A7S='R' @ A8S='I' 
A9S=""ROOT#',DD,': ',A,'=',MZ.6DE" 

DEF FNF2(U1,V1) =SQRT (U1*U1+V1*V1) 

DISP "POLYNOMIAL ROOT FINDER" @ WAIT .5 


INPUT ORDER, COEFFICIENTS, TOLERANCES **** 
AND MAXIMUM NUMBER OF ITERATIONS SRS*X 


"ORD': 
INPUT "ORDER OF POLYNOMIAL? "eN 

IF N<=1l OR N#IP(N) THEN DISP “INVALID ORDER" @ GOTO 'ORD' 
ON ERROR GOTO 'MEM' 

DIM A@(N) ,R@(N,2) ,L(2) ,C9(N) 

FOR I=N TO @ STEP -1 

DISP ‘A(';STRS(I);') ="; @ INPUT A@(I) 
NEXT I 

'TOL': 

El=.08900000001 @ SCI @ 

INPUT ‘TOL. FOR ROOTS=',STRS(E1) ;El 
INPUT 'TOL. FOR FCN=',STRS (E1*10@) ;E2 
STD 

‘IT's: INPUT "MAX # OF ITERATIONS=' ; 19 
IF IP(I@)#I@ OR NOT I@ THEN 'IT' 

DISP "LOOKING FOR ROOTS...' @ WAIT .5 
J=N @ K-G @x=¢ @ y=¢ 

'FINDR': 


IS ZERO A ROOT? **** 
IF NOT A@(@) THEN ‘FAR‘ 
INPUT NEW GUESS IF ITERATION LIMIT EXCEEDED **** 


‘LOOP': K=K+l 
IF K>I®@ THEN DISP 'NO CONVERGENCE' @ WAIT 1 @ GOSUB 'NUROOT! 


CALCULATE P, P', AND P'' AT Z(x,y) **** 


R=X*X+Y*Y @ D=X+X 

D@=G @ D1=8 @ CH=G @ C1=9G 
B1=A@(J) @ BO=AG(J-1)+D*AG (J) 
FOR I=J-2 TO @ STEP -1 

IF I=@ THEN D=X 

IF I>J-4 THEN 38@ 

V=D1*R @ D1=DGd 
DG=C1+D*DG-V 

V=C1*R @ C1=CG 
C8=B1+D*C8-V 

V=B1*R @ B1=BG 
BO=A9 (I) +D*BO-V 

NEXT I 
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kkk 


8430 
8448 
8450 
84608 
0478 
G48 
G49 
6580 
0519 
05208 
G53@ 
G54 
8558 


kakKK 
kkkk 


8568 
0578 


kkk 


8586 
8598 


kKKK 


8690 
8618 
8626 


kkk 
kkk 


8639 
6649 
8659 
8660 
8676 
8680 


kkKK 


8698 
0700 
8719 
0728 
0736 
0746 
8758 
0768 


LISTING 


BEGIN CALCULATION OF LAGUERRE STEP **** 


P(1,1)=BO @ P(1,2)=B1*Y 

P(2,1)=B1-2*Y*Y*Cl @ P(2,2)=2*Y*CO 
P(3,1)=2*CO-8*Y*Y*DG @ P(3,2) =2*Y* (3*C1-4*Y*Y*D1) 
CALL MULT(P(1,1),P(1,2),P(3,1) ,P (3,2) ,S1,S2) 
$1l=-S1*J*(J-1) @ $2=-S2*J* (J-1) 

CALL MULT (P(2,1),P(2,2) ,P (2,1) ,P (2,2) ,S3,84) 
$3=S3*(J-1)*2 @ S4=S4*(J-1)%*2 

CALL ADD(S3,84,81,S2,S1,S82) 

CALL SQAROOT(S1,S2,S1,S82) 

CALL ADD(P(2,1),P(2,2),S1,S2,L1,L2) 

CALL ADD(P(2,1),P(2,2) ,-Sl1,-S2,L3,L4) 

CALL DIVID(P(1,1),P(1,2) ,L1,L2,L1,L2) 

CALL DIVID(P(1,1),P(1,2) ,L3,L4,3,L4) 


CHOOSE SIGN OF DENOMINATOR TO **** 
PRODUCE SMALLER LAGUERRE STEP **** 


IF FNF2(L1,L2)<FNF2(L3,L4) THEN S1=L1 @ S2=L2 ELSE S1=L3 @ S2=L4 
L(1)=-3*S1 @ L(2)=-J*S2 


FORCE REAL ROOT **** 


IF ABS(L(2))<.@@@1*FNF2(L(1),L(2)) THEN L(2)=@8 @ Y=@ 
X=X+L(1) @ Y=Y¥+L(2) 


CHECK FOR VALID ROOT **** 


IF FNF2(L(1),L(2))>ABS(E1l) THEN ‘LOOP’ 
IF FNF2(P(1,1),P(1,2))<ABS(E2) THEN 'FAR' 
DISP 'INVALID ROOT FOUND' @ WAIT 1 @ GOSUB 'NUROOT' @ GOTO ‘LOOP’ 


FOUND VALID ROOT(S) - IF COMPLEX **** 
ROOT ASSUME CONJUGATE IS A ROOT eAAS 


‘FAR': 

RO(J,1)=X @ RO(J,2)=¥ @ IF Y THEN GOSUB 'ROOTI' 

DISP '# OF ROOTS FOUND =';N-J+1 

IF NOT FNF2(R@(J,1),R@(J,2)) THEN GOSUB 'ROOT@' ELSE GOSUB 'DEFLATE' 
J=J-1 @ K=@ @ IF NOT J THEN 'DR' 

IF J=l1 THEN R@(1,1)=-AG@(G)/AG(1) ELSE K=@ @ X=@ @ ‘Ys @ @ GOTO 'FINDR' 


DISPLAY ROOTS **** 


DISP '# OF ROOTS FOUND =';N 

"DR': 

FOR I=N TO 1 STEP -1 

DISP USING A9S;N-I+1,A7$,R0(I,1) @ GOSUB 'WAIT' 
DISP USING A9S;N-1+1,A8S,RG (1,2) @ GOSUB ‘'WAIT' 


NEXT I 
"DONE': PUT "#43" @ END 
'WAIT': IF KEYS='' THEN ‘WAIT’ ELSE RETURN 
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kkk 
kkKK 
ekKK 


8776 
8786 
87998 
8890 
8818 
8828 
88306 
8846 
0856 
8869 
8870 
8880 
8899 
8988 
8918 


kekk 


8928 
G93 
8940 


xkkkk 
xkxkkk 


8950 
6960 
0978 


kekk 


0986 
8990 
1688 
1919 
1829 
1938 
1949 
1858 
1968 
1876 
10898 


kkk 


1999 
1169 
1114 


LISTING 


DEFLATION ROUTINES - IF ROOT IS REAL, DEFLATE BY **** 
LINEAR FACTOR (x-r). IF ROOT IS COMPLEX, DEFLATE **** 
BY BINOMIAL X*2-2Re(r)X+!iri!*2. Keke 


"DEFLATE': 

IF Y THEN 'DEFLATEIL' 

C9 (J-1) =AG (J) 

FOR I=J-1 TO 1 STEP -1 

C9 (I-1) =AG(1I)+C9 (I) *RO(J,1) 

NEXT I 

FOR I=@ TO J-1 @ A®(I)=C9(I) @ NEXT I 

RETURN 

'DEFLATEL': 

C9(J-1)=AG(J+1) @ C9(J)=8 @ IF J=1 THEN RETURN 
FOR I=J-1 TO 1 STEP -1l 

C9 (I-1) =AG(I+1)+C9 (I) *RO(J,1)*2-(RO(T,1) *2+RO(J,2) ~2) *C9 (I+1) 
NEXT I 

FOR I=@ TO J-1 @ AG(I)=C9(I) @ NEXT I 

RETURN 


DEFLATE FOR ZERO ROOT **** 


‘ROOTO': 
FOR I=@ TO J-1 @ AG(I)=AG(I+1) @ NEXT I 
RETURN 

SET NEXT ROOT TO COMPLEX CONJUGATE **** 
OF COMPLEX ROOT JUST FOUND zeee 
*ROOTI': 
J=J-1 @ RO(J,1)=X @ RO(J,2)=-Y 

RETURN 


ASK FOR NEW GUESS IF FAILURE TO CONVERGE **** 


"NUROOT' : 

DISP 'NEW GUESS? (Y/N)' 

'WO': AlS=KEYS @ IF NOT POS(A@S,UPRCS(A1$)) THEN 'wo' 
IF UPRCS(A1S)='N' THEN 'DONE' 

INPUT 'NEW u =',STRS(X);X 

INPUT 'NEW v =',STRS(Y) ;Y 

K=@ 

RETURN 

*MEM': 

IF ERRN#24 THEN DISP ERRMS @ GOTO 'DONE' 

DISP "LOW MEM - REDUCE ORDER' @ WAIT 1 @ GOTO 'ORD' 


ADDITION OF COMPLEX NUMBERS **** 
SUB ADD(U1,V1,U2,V2,U,V) 


U=U1+U2 @ V=V1lt+v2 
END SUB 
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kkkk 


1129 
1138 
1148 


kkKK 


115@ 
11698 
1179 
118@ 
1198 
1290 


kkk 


1219 
1226 
12398 
1246 
1258 
12698 
1276 
1288 


LISTING 


MULTIPLICATION OF COMPLEX NUMBERS **** 
SUB MULT(U1,V1,U2,V2,U,V) 

U=U1*U2-V1*V2 @ V=U1*V2+U2*V1 

END SUB 

DIVISION OF COMPLEX NUMBERS **** 

SUB DIVID(U1,V1,U2,V2,U,V) 
D5=U2*U2+V2*V2 @ IF NOT D5 THEN U=@ @ V=@8 @ GOTO 12989 
U=Z1/D5 
V=Z2/D5 

END SUB 

SQUARE ROOT OF A COMPLEX NUMBER **** 
SUB SQAROOT(U1,V1,U,V) 

A2=SQR ((SQR(U1*U1+V1* Vf +ABS (U1) ) /2) 

IF NOT A2 THEN U=9 @ V=8 @ GOTO 1286 
B2=V1/(2*A2) 

IF Ul>=@ THEN U=A2 @ V=B2 @ GOTO 1286 
U=ABS (B2) 

IF B2>=@ THEN V=A2 ELSE V=-A2 

END SUB 
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(AI mackano 


00071-90064 Printed in USA 


ERRATA: HP-71 MATH SOLUTIONS BOOK 


Page 57: : 
5 big 
3) Enter the coefficients, 
starting with the highest 
order term. 


\4 Calculation of roots begins. 


\ 8) Display real and imaginary 
parts of each root. 


Page 58: 


0253 J=N @ K=O @ X=9 @ ysg 
@270 IF NOT AG(G) THEN 'FAR' 


.Page 59: 


5 [ENDLINE] 
-45 [ENDLINE] 
225 [ENDLINE] 

-425 [ENDLINE}) 


176 [ENDLINE] 


370 (ENDLINE]} 
-599 [ENDLINE] 


<any key> 
<any key> 
<any key> 
<any key> 
<any key> 


A(5)=? 
A(4)=? 
A(3)=? 
A(2) =? 
A(1)=? 
A(@) =? 


5/16/84 


TOL. FOR ROOTS=1.E-1¢6 


# OF ROOTS FOUND= 
# OF ROOTS FOUND= 
# OF ROOTS FOUND= 
# OF ROOTS FOUND= 


ROOT# 
ROOT # 
ROOT# 
ROOT # 
ROOT# 


- ROOT# 


WWND RH 


oe e808 60 @¢ ce o¢ 


6728 DISP USING A9$;N-I+1,A7S,R0(1I,1) @ GOSUB 'WAIT' 
@73@ DISP USING A9S N-I+1,A8$,R0(1I,2) @ GOSUB ‘WAIT’ 
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1220 A2=SQR((SQR(U1*U1+V1*V1) +ABS (U1) )/2) 


OP wWnh 


R= 1.000300E+009 
T= 1.000300E+006 
R= 1.009000E+900 
I=-1.003090E+000 


. R=-1. AO9SG0E+O60G 
‘T= 0. 08900GE+000 


ERRATA: HP-71 MATH SOLUTIONS BOOK 
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\ g200 IF NOT FLAG(1) EXOR FLAG(2) THEN DISP 


“0518 IF FLAG (2) THEN CALL C2R(X@,YB) 
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\ 4) Key in the number of points. 9 [ENDLINE] 


Xe) Key in the interval length 


for the partitions. 
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~ 9580 R=A(G)+A(N) 


Page 27: 
\@99@ MS=KEYS @ IF MS="" 
\ G16 AS=KEYS @ IF AS="" 
'\ 6238 AS=KEYS @ IF AS="" 


Page 28:: 
9418 AS=KEY$ @ IF AS="" 


Page 34: 
° \3) Key in the number of 


Page 36: 


\N3) The order is 2. 
Ng) Choose solve option. 
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- 25 [ENDLINE] 


THEN 99 
THEN 168 
THEN 239 ELSE "F" 


THEN 416 


rows. 5 [ENDLINE] 


2 
[ENDLINE] 


S 


c 6) Frequency domain output. 


<any key> 
<any key> 


5/16/84 


"I/O INCORRECT’ @ WAIT 2 ELSE ‘Al 


Interval length=@ 


E (8) = 


A(1,1)= 6.0906 


R(1,1)= 0.9606 


BR(1)= 4.0000 


DC TERM=@ 
MAX FREQ. =@ 

FREQ DOMAIN OUTPUT: 
1R= 1. 00001GE+906 


